Introducing ordination: distance-based methods

Gabriel Singer

24.02.2022

A short intro to classical ordination methods - in R ;-)

Ecological datasets have one or both of these datasets (if not more):

A short intro to classical ordination methods - in R ;-)

The classical ordination methods always target dimension reduction, and

Note that ordination is a word almost exclusively used by ecologists. Outside most of these methods are known as scaling techniques. The phrase gradient analysis makes a bit more sense in ecology ;-)

Specifically PCA and RDA are quite useful outside of the classical species or species~environment framework of ecology.

Principal component analysis (PCA)

PCA can be understood as a MLR on theoretical (latent) instead of observed dependent variables.

The regression coefficients in PCA are factor loadings.

The predicted values of the theoretical dependent are scores.

For instance, in a PCA on this matrix of n objects times p variables PC1 is computed as:

\[ PC_1=f_1X_1+f_2X_2+...+f_pX_p \] There must be as many PCs as original variables. The factor loading matrix translates p \(X\)-variables into p \(PCs\).

PCA can be understood as a rotation of the original coordinate system made up of p \(X\)-variables into a new coordinate system defined by p \(PCs\).


Principal component analysis (PCA) in R

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
library(shape) # nice arrows
mara<-read.table(file="data/MaraRiver.txt",header=TRUE) # water chemistry in 54 streams, 3 types of landuse

wc<-log(mara[-c(5,6),9:26]) # delete two cases with NA data, log-transform concentration data
landuse<-mara$landuse[-c(5,6)]
plot(wc) # check potential for PCA by correlation plot


pca<-prcomp(scale(wc),retx=T,center=F,scale.=F)
pca<-prcomp(wc,retx=T,center=T,scale.=T)   # equivalent to line above

summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.7152 1.9490 1.22235 1.11919 0.88947 0.84667 0.80815
## Proportion of Variance 0.4096 0.2110 0.08301 0.06959 0.04395 0.03982 0.03628
## Cumulative Proportion  0.4096 0.6206 0.70362 0.77321 0.81716 0.85699 0.89327
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.69942 0.68698 0.57278 0.44077 0.34335 0.32238 0.29800
## Proportion of Variance 0.02718 0.02622 0.01823 0.01079 0.00655 0.00577 0.00493
## Cumulative Proportion  0.92045 0.94667 0.96489 0.97569 0.98224 0.98801 0.99294
##                           PC15    PC16    PC17    PC18
## Standard deviation     0.26124 0.20738 0.09562 0.08134
## Proportion of Variance 0.00379 0.00239 0.00051 0.00037
## Cumulative Proportion  0.99674 0.99912 0.99963 1.00000

pca$sdev # stdevs of PCs (squares are eigenvalues)
##  [1] 2.71519185 1.94904622 1.22234767 1.11918666 0.88947376 0.84666850
##  [7] 0.80815199 0.69941920 0.68698308 0.57277816 0.44076721 0.34334955
## [13] 0.32237891 0.29800143 0.26124205 0.20738029 0.09561853 0.08134109

head(scores<-pca$x) # site scores on all PCs
##         PC1      PC2        PC3        PC4        PC5         PC6        PC7
## 1 -5.170835 4.016387 -0.5991588  1.2909046 -0.1010020  0.16667942  0.6343563
## 2 -4.109584 3.447343  0.1318830 -2.3348590 -0.7292225 -0.47398838  0.3727337
## 3 -3.448190 4.451837 -0.9213255  0.5202531 -0.6210845  0.55669779 -1.1869106
## 4 -4.026131 2.619122  0.3155890  0.2170746  0.2018600 -0.04762554 -0.7586118
## 7  2.232793 1.480272  0.4206951  0.6272105  1.5808947 -0.34864745 -0.5814405
## 8  5.913178 1.424879 -3.2684397 -2.9349723 -1.0034926  1.03815652 -0.2327136
##          PC8        PC9       PC10       PC11        PC12        PC13
## 1 -0.2298109 -0.9806756 -0.3871041  0.3571909 -0.04236249 -0.03051565
## 2 -0.1099423 -1.3833508 -1.1315837 -0.2817962 -0.63408704  0.33104112
## 3 -0.6201881  0.1749492 -0.5086871  0.6001108 -0.52214446  0.20520087
## 4 -1.0545257  0.6203206 -0.6279348  0.1036275  0.48580457 -0.21386990
## 7 -0.9861408  1.3964450 -0.5847160 -1.2657650 -0.48722114 -0.36661223
## 8 -0.2793900  0.5258090  0.7129195 -0.1997592 -0.26750430 -0.60819198
##         PC14        PC15        PC16        PC17        PC18
## 1 -0.4348322 -0.06170742  0.03812678 -0.06946439  0.04456263
## 2  0.1136519  0.13224238  0.15079191 -0.27204153 -0.18865812
## 3  0.5252942 -0.62886628 -0.05060063  0.10903471  0.06976005
## 4  0.1518358  0.38178030 -0.12215586 -0.07002405  0.07068122
## 7 -0.3610469  0.48773223  0.09650913  0.06186595  0.03892577
## 8 -0.5882226 -0.29698046 -0.07645784 -0.11833869 -0.12134436

head(scale(wc) %*% pca$rotation) # to manually compute scores from variables and loadings
##         PC1      PC2        PC3        PC4        PC5         PC6        PC7
## 1 -5.170835 4.016387 -0.5991588  1.2909046 -0.1010020  0.16667942  0.6343563
## 2 -4.109584 3.447343  0.1318830 -2.3348590 -0.7292225 -0.47398838  0.3727337
## 3 -3.448190 4.451837 -0.9213255  0.5202531 -0.6210845  0.55669779 -1.1869106
## 4 -4.026131 2.619122  0.3155890  0.2170746  0.2018600 -0.04762554 -0.7586118
## 7  2.232793 1.480272  0.4206951  0.6272105  1.5808947 -0.34864745 -0.5814405
## 8  5.913178 1.424879 -3.2684397 -2.9349723 -1.0034926  1.03815652 -0.2327136
##          PC8        PC9       PC10       PC11        PC12        PC13
## 1 -0.2298109 -0.9806756 -0.3871041  0.3571909 -0.04236249 -0.03051565
## 2 -0.1099423 -1.3833508 -1.1315837 -0.2817962 -0.63408704  0.33104112
## 3 -0.6201881  0.1749492 -0.5086871  0.6001108 -0.52214446  0.20520087
## 4 -1.0545257  0.6203206 -0.6279348  0.1036275  0.48580457 -0.21386990
## 7 -0.9861408  1.3964450 -0.5847160 -1.2657650 -0.48722114 -0.36661223
## 8 -0.2793900  0.5258090  0.7129195 -0.1997592 -0.26750430 -0.60819198
##         PC14        PC15        PC16        PC17        PC18
## 1 -0.4348322 -0.06170742  0.03812678 -0.06946439  0.04456263
## 2  0.1136519  0.13224238  0.15079191 -0.27204153 -0.18865812
## 3  0.5252942 -0.62886628 -0.05060063  0.10903471  0.06976005
## 4  0.1518358  0.38178030 -0.12215586 -0.07002405  0.07068122
## 7 -0.3610469  0.48773223  0.09650913  0.06186595  0.03892577
## 8 -0.5882226 -0.29698046 -0.07645784 -0.11833869 -0.12134436
head(loadings<-pca$rotation) # the variable loadings
##              PC1        PC2        PC3        PC4         PC5         PC6
## cond -0.25722323  0.1457565 -0.3356613 -0.1583734 -0.04126487 -0.40451583
## pH    0.05614221  0.1125082  0.5222340 -0.3741856 -0.17110129 -0.25878527
## turb -0.13203932 -0.3855016 -0.1073370 -0.3716724 -0.22025887 -0.01157911
## Alk  -0.13633672  0.3519383 -0.2457525  0.0313204  0.09063681 -0.44809661
## TSS  -0.11636881 -0.3402503 -0.1025929 -0.4995496 -0.25110405  0.04368404
## PO4  -0.08826034 -0.2141495  0.4323275 -0.1977583  0.60913902 -0.12539930
##              PC7          PC8         PC9         PC10        PC11        PC12
## cond  0.10974673 -0.050428262  0.17413681  0.013940888  0.64808704 -0.08137442
## pH   -0.58257813  0.096086825  0.26386503 -0.205980904  0.05773032 -0.05873929
## turb  0.08448051 -0.011737527 -0.06980739  0.001650434 -0.31037110 -0.49543374
## Alk   0.04261110 -0.306722599  0.20931630 -0.141917638 -0.60093241 -0.08602031
## TSS   0.13099111 -0.073853481 -0.15289911 -0.157532776 -0.05344048  0.39788229
## PO4   0.40877114 -0.002124487  0.17867170 -0.021848936 -0.06745641  0.06817771
##             PC13        PC14        PC15          PC16         PC17
## cond  0.07842483 -0.11147245  0.19970689 -0.2485760664 -0.059467443
## pH   -0.05666025 -0.03676865 -0.04291062  0.0310385860 -0.004463366
## turb  0.38994704 -0.19755992 -0.16242245 -0.2347958749  0.090676524
## Alk  -0.18762446  0.11999960  0.01473726  0.0008990899 -0.068827445
## TSS  -0.36793380  0.35906128  0.20457551 -0.0570037711 -0.064652896
## PO4   0.03969603 -0.21939294  0.26938989  0.0462139650 -0.019837594
##              PC18
## cond  0.118215868
## pH    0.018664958
## turb  0.029801918
## Alk   0.065155143
## TSS   0.009829404
## PO4  -0.007232622

# how many axes should be kept?)
screeplot(pca,npcs=length(pca$sdev),type = "lines") # plots eigenvalues vs. component #
abline(h=1,col="red")

# -> the first 3-4 PCs seem useful, and just PC1 and PC2 alone are already explaining a lot of overall variance
# eigenvalue of PC5<1, so PC5 contributes less than one original variable (Kaiser-Guttman criterion)

# for a DISTANCE BIPLOT (focus is on sites, "scaling 1")
# each principal component has variance given by eigenvalue, loadings remain unscaled
plot(scores[,1:2],asp=1,pch=21,bg=landuse) # note asp=1
arrows<-loadings*7 # with extension factor
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="darkgreen")
text(x=arrows[,1]*1.3,y=arrows[,2]*1.2,labels=names(wc),cex=0.7)


biplot(pca,scale=0)

# in this plot:
# 1) distances among sites are approximating true Euclidean distances in multivariate space
# 2) angles between arrows do not reflect correlations among variables
# 3) projecting site on descriptor at right angle gives its appr. descriptor value

# for a CORRELATION BIPLOT (focus is on variables, "scaling 2")
# each principal component is weighted by 1/sqrt(eigenvalue), so it has variance 1
var(scores[,1]/pca$sdev[1]) # just demo
## [1] 1
plot(scores[,1]/pca$sdev[1],scores[,2]/pca$sdev[2],pch=21,bg=landuse,asp=1)
# loadings are weighted by sqrt(eigenvalues)
arrows<-loadings*matrix(pca$sdev,nrow=nrow(loadings),ncol=ncol(loadings),byrow=TRUE)
arrows<-arrows*2 # choose extension factor
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="purple")
# as alternative just compute correlation of scores with original data ("structure coefficients")
(structure<-cor(wc,scores))
##             PC1         PC2         PC3          PC4         PC5          PC6
## cond -0.6984104  0.28408615 -0.41029477 -0.177249446 -0.03670402 -0.342490809
## pH    0.1524369  0.21928376  0.63835148 -0.418783552 -0.15219011 -0.219105339
## turb -0.3585121 -0.75136050 -0.13120308 -0.415970767 -0.19591449 -0.009803667
## Alk  -0.3701803  0.68594395 -0.30039498  0.035053376  0.08061906 -0.379389286
## TSS  -0.3159637 -0.66316347 -0.12540419 -0.559089276 -0.22335047  0.036985903
## PO4  -0.2396438 -0.41738723  0.52845449 -0.221328417  0.54181317 -0.106171637
## NH4  -0.4654270 -0.02194250 -0.53020009 -0.186343168  0.45612324  0.162850758
## NO2  -0.5188878 -0.52230884 -0.01574993  0.005211432  0.23769161 -0.002141302
## NO3  -0.5730405 -0.65160683  0.03213869  0.382790600 -0.08425720  0.030755129
## Br   -0.9313245 -0.07779131 -0.01608737  0.019147954 -0.14787110  0.071294662
## Cl   -0.8345065 -0.10268275  0.30020535  0.279338280  0.10033310 -0.071715149
## Fl   -0.8480805 -0.07571208  0.09535793 -0.020000354  0.06993200 -0.351671172
## SO4  -0.7986422  0.21769902  0.23421111  0.278802688 -0.11423549  0.046932670
## Na   -0.8811191  0.16618532  0.01657297  0.038583410 -0.22572697 -0.089087312
## K    -0.7557029  0.41508974  0.24253708 -0.075909350 -0.00631684  0.223742281
## Ca   -0.6701345  0.61799716 -0.01657030 -0.233382943  0.03311222  0.221873328
## Mg   -0.7070641  0.54017541  0.09521808 -0.211787735  0.00985287  0.346221391
## TDN  -0.6378635 -0.61620157 -0.10341460  0.256075454 -0.06798394  0.064818416
##              PC7          PC8          PC9          PC10         PC11
## cond  0.08869204 -0.035270495  0.119629039  0.0079850361  0.285655516
## pH   -0.47081167  0.067204970  0.181270813 -0.1179813623  0.025445632
## turb  0.06827309 -0.008209452 -0.047956499  0.0009453327 -0.136801402
## Alk   0.03443624 -0.214527675  0.143796760 -0.0812873231 -0.264871300
## TSS   0.10586072 -0.051654543 -0.105039104 -0.0902313332 -0.023554813
## PO4   0.33034921 -0.001485907  0.122744436 -0.0125145935 -0.029732574
## NH4  -0.32682252  0.325631695 -0.037652130 -0.1272200457 -0.032943687
## NO2  -0.36795338 -0.451185072 -0.141524231  0.1951870779  0.047492681
## NO3  -0.05158254  0.000133211  0.277017753 -0.0433056814 -0.016570634
## Br    0.14345062 -0.019023291 -0.004000929 -0.1725441479  0.066895666
## Cl    0.01996508  0.024214082 -0.120525858 -0.1886080680  0.048860822
## Fl    0.03972002  0.187589873 -0.200265391  0.1477822419 -0.005407743
## SO4  -0.07531681 -0.041088827 -0.301983239 -0.1831854495 -0.043669377
## Na   -0.08523333  0.179791964 -0.023950118  0.2200084209 -0.091597546
## K     0.07557981  0.121199537  0.072790327  0.2593618792 -0.007213583
## Ca    0.04548095 -0.165843239  0.100610952 -0.0774684592  0.016715739
## Mg    0.01839427 -0.106342119  0.063503614  0.0234741312 -0.013933336
## TDN  -0.09445248  0.037539759  0.334130611  0.0077866617 -0.019725516
##             PC12         PC13          PC14         PC15          PC16
## cond -0.02793987  0.025282511 -0.0332189485  0.052171838 -0.0515497755
## pH   -0.02016811 -0.018266069 -0.0109571105 -0.011210059  0.0064367908
## turb -0.17010695  0.125710701 -0.0588731392 -0.042431573 -0.0486920355
## Alk  -0.02953504 -0.060486168  0.0357600539  0.003849992  0.0001864535
## TSS   0.13661271 -0.118614098  0.1070007750  0.053443725 -0.0118214583
## PO4   0.02340879  0.012797165 -0.0653794093  0.070375967  0.0095838652
## NH4  -0.01376245 -0.020251926 -0.0060367305  0.016144920  0.0028392903
## NO2  -0.02645191 -0.020466386  0.0025810284  0.009774220  0.0245358000
## NO3   0.06814171 -0.023152938 -0.0062065455 -0.018023245 -0.0252696403
## Br   -0.10337675 -0.088756432 -0.0618928002 -0.030769421  0.1430478324
## Cl   -0.07884890  0.093133026  0.2026581905  0.002730477 -0.0151599829
## Fl    0.10531097 -0.019250758 -0.0169508010 -0.162688990  0.0023135739
## SO4   0.03177807 -0.044383121 -0.1385124348  0.063800533 -0.0813112489
## Na    0.05382217  0.106875862  0.0049270489  0.137864470  0.0702378702
## K    -0.13408408 -0.158024615  0.0489002643  0.007774656 -0.0586163075
## Ca    0.08575336  0.066801075 -0.0031824718 -0.045231364 -0.0142660407
## Mg    0.06245134  0.090206616  0.0004852794 -0.046904944  0.0073118804
## TDN   0.03187638 -0.003105378 -0.0152069323 -0.017003161 -0.0228256456
##               PC17          PC18
## cond -0.0056861896  0.0096158080
## pH   -0.0004267806  0.0015182281
## turb  0.0086703562  0.0024241206
## Alk  -0.0065811793  0.0052997906
## TSS  -0.0061820150  0.0007995345
## PO4  -0.0018968416 -0.0005883094
## NH4   0.0050891691  0.0049022516
## NO2   0.0027370353 -0.0011256331
## NO3   0.0441001155  0.0418426505
## Br    0.0009045142 -0.0009829097
## Cl   -0.0034330802 -0.0031772767
## Fl   -0.0013409888 -0.0022312652
## SO4  -0.0082417516 -0.0023487676
## Na    0.0102682671 -0.0035812642
## K     0.0081296716 -0.0014049273
## Ca    0.0495490253 -0.0377676853
## Mg   -0.0458882844  0.0395256212
## TDN  -0.0464454202 -0.0410550878
structure<-2*structure
Arrows(x0=0,y0=0,x1=structure[,1],y1=structure[,2],col="red")
text(x=arrows[,1]*1.3,y=arrows[,2]*1.2,labels=names(wc),cex=0.7)


biplot(pca,scale=1)

# in this plot
# 1) distances among sites are not approximating true Euclidean distances in multivariate space
# 2) angles between arrows reflect correlations among variables (NOT proximity of arrow heads)
# 3) projecting site on descriptor at right angle gives its appr. descriptor value

# PCA using the rda function from the vegan package
pca2<-rda(X=wc,scale=TRUE)
summary(pca2,scaling=1)
## 
## Call:
## rda(X = wc, scale = TRUE) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total              18          1
## Unconstrained      18          1
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Eigenvalue            7.3723 3.7988 1.49413 1.25258 0.79116 0.71685 0.65311
## Proportion Explained  0.4096 0.2110 0.08301 0.06959 0.04395 0.03982 0.03628
## Cumulative Proportion 0.4096 0.6206 0.70362 0.77321 0.81716 0.85699 0.89327
##                           PC8     PC9    PC10    PC11     PC12     PC13
## Eigenvalue            0.48919 0.47195 0.32807 0.19428 0.117889 0.103928
## Proportion Explained  0.02718 0.02622 0.01823 0.01079 0.006549 0.005774
## Cumulative Proportion 0.92045 0.94667 0.96489 0.97569 0.982236 0.988010
##                           PC14     PC15     PC16      PC17      PC18
## Eigenvalue            0.088805 0.068247 0.043007 0.0091429 0.0066164
## Proportion Explained  0.004934 0.003792 0.002389 0.0005079 0.0003676
## Cumulative Proportion 0.992944 0.996735 0.999124 0.9996324 1.0000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.504409 
## 
## 
## Species scores
## 
##          PC1      PC2      PC3      PC4      PC5      PC6
## cond -1.4159  0.80230 -1.84762 -0.87175 -0.22714 -2.22662
## pH    0.3090  0.61929  2.87459 -2.05967 -0.94181 -1.42446
## turb -0.7268 -2.12196 -0.59083 -2.04584 -1.21239 -0.06374
## Alk  -0.7505  1.93721 -1.35272  0.17240  0.49890 -2.46651
## TSS  -0.6405 -1.87288 -0.56471 -2.74973 -1.38218  0.24045
## PO4  -0.4858 -1.17877  2.37971 -1.08854  3.35295 -0.69025
## NH4  -0.9435 -0.06197 -2.38757 -0.91648  2.82267  1.05873
## NO2  -1.0519 -1.47508 -0.07092  0.02563  1.47093 -0.01392
## NO3  -1.1617 -1.84024  0.14473  1.88265 -0.52142  0.19995
## Br   -1.8880 -0.21969 -0.07244  0.09417 -0.91508  0.46350
## Cl   -1.6918 -0.28999  1.35187  1.37385  0.62090 -0.46624
## Fl   -1.7193 -0.21382  0.42941 -0.09837  0.43277 -2.28630
## SO4  -1.6191  0.61482  1.05469  1.37121 -0.70693  0.30512
## Na   -1.7863  0.46933  0.07463  0.18976 -1.39689 -0.57918
## K    -1.5320  1.17228  1.09218 -0.37334 -0.03909  1.45461
## Ca   -1.3585  1.74532 -0.07462 -1.14783  0.20491  1.44246
## Mg   -1.4334  1.52554  0.42878 -1.04162  0.06097  2.25087
## TDN  -1.2931 -1.74025 -0.46569  1.25944 -0.42071  0.42140
## 
## 
## Site scores (weighted sums of species scores)
## 
##         PC1      PC2      PC3       PC4      PC5       PC6
## 1  -0.93940  0.72967 -0.10885  0.234522 -0.01835  0.030281
## 2  -0.74660  0.62629  0.02396 -0.424180 -0.13248 -0.086111
## 3  -0.62644  0.80878 -0.16738  0.094516 -0.11283  0.101137
## 4  -0.73144  0.47582  0.05733  0.039436  0.03667 -0.008652
## 7   0.40564  0.26892  0.07643  0.113947  0.28721 -0.063340
## 8   1.07426  0.25886 -0.59379 -0.533204 -0.18231  0.188605
## 9  -0.32892  0.25033 -0.45377  0.617791 -0.24044 -0.122356
## 10 -0.57184  0.21769  0.03691 -0.379616  0.08341 -0.294496
## 11  0.13734  0.33630 -0.26902  0.054280  0.12820 -0.050003
## 12  0.78569  0.17837  0.01844  0.115072 -0.02183  0.091719
## 13 -0.22093 -0.36549 -0.16998 -0.099846  0.09512 -0.143509
## 14  1.30998 -0.67125 -0.57884 -0.023233  0.14289 -0.479781
## 15 -0.43427 -0.59259 -0.01832 -0.068894 -0.07069  0.077409
## 16 -0.74333 -0.01680 -0.19101  0.152920  0.30242 -0.080449
## 17 -0.84366 -0.21000 -0.04445 -0.008928  0.07997  0.014100
## 18 -0.48203 -0.01707 -0.11851 -0.044274 -0.16607 -0.277649
## 19 -0.28152 -0.26355  0.12280 -0.173659 -0.05794  0.076838
## 20  0.01911 -0.20131 -0.17851  0.014238 -0.28113  0.127905
## 21 -0.40327 -0.31194  0.14452 -0.048581 -0.05699  0.013663
## 22 -0.18831 -0.29426  0.20016 -0.001522 -0.30260  0.030065
## 23 -0.42134 -0.16552  0.05614 -0.096723  0.05274  0.079374
## 24 -0.67474 -0.45065 -0.01582 -0.035151  0.15207  0.097047
## 25 -0.07187 -0.38842  0.37212 -0.278347 -0.02696 -0.088559
## 26  0.16221 -0.46968  0.08605 -0.082842  0.04253  0.100862
## 27  0.02137 -0.26587 -0.13515  0.126573 -0.12436  0.198872
## 28  0.04144 -0.47296 -0.01313 -0.042046  0.11692  0.111934
## 29  0.04714 -0.40758  0.03664  0.034209  0.05919  0.114249
## 30 -0.18686 -0.40929 -0.27928  0.050715 -0.07947  0.204339
## 31  0.37961  0.30976  0.11739  0.216706 -0.20004 -0.433454
## 32  0.17514  0.28779 -0.02441  0.098778  0.05891  0.112364
## 33  0.05528 -0.03829  0.12601  0.556220  0.21520 -0.184383
## 34 -0.12436  0.23716  0.22438 -0.342063  0.07600 -0.046932
## 35 -0.05549 -0.18388  0.16587 -0.045118  0.03482 -0.030455
## 36  0.20401  0.36785 -0.01204 -0.020236 -0.01134  0.003642
## 37  0.54387  0.14867 -0.04013  0.095238 -0.18598  0.179905
## 38  0.35917  0.14769  0.17737 -0.004264  0.22119  0.093820
## 39  0.42440  0.35483  0.14198  0.045518  0.09675  0.190421
## 40  0.33889  0.19385  0.14972 -0.003335  0.25193  0.134582
## 41  0.71629  0.08022  0.68924  0.002355 -0.12205 -0.199364
## 42  0.78687  0.25324  0.02979  0.145646 -0.20004  0.072693
## 43  0.22158  0.29242  0.32534 -0.093631  0.06309 -0.015564
## 44  0.61384  0.01066  0.13364  0.126831  0.28036  0.112728
## 45  0.43087  0.68355 -0.25570 -0.420474  0.16347 -0.019426
## 46 -0.34968 -0.38830 -0.17466 -0.086047  0.04080 -0.195030
## 47 -0.12295 -0.09128 -0.05881 -0.061097  0.09589 -0.145553
## 48  0.21811  0.18380  0.03419  0.164258  0.12432  0.150992
## 49 -0.11584 -0.02878 -0.15463 -0.137653 -0.15426  0.067831
## 50 -0.16384 -0.28048 -0.04759  0.027258 -0.12822  0.077030
## 51 -0.28502 -0.08784  0.14403  0.097088 -0.02838 -0.002284
## 52  0.43786 -0.05135  0.34622  0.043325 -0.44642 -0.131500
## 53  0.11174 -0.24041 -0.02979  0.099997  0.10370  0.125014
## 54  0.09225 -0.33766  0.09691  0.187527 -0.05459  0.119429
scores(wc,scaling=0)
##        cond       pH       turb      Alk      TSS       PO4        NH4
## 1  5.434159 2.008214 -3.9120230 7.677864 2.529802 -14.47716 -10.688750
## 2  5.324959 2.083185  3.5695327 7.677864 4.085042 -13.99943 -10.654264
## 3  5.326419 2.064328 -3.9120230 7.677864 2.358155 -15.80328 -10.490475
## 4  5.202907 2.074429 -3.9120230 7.677864 3.294588 -13.81551 -10.785377
## 7  4.169761 2.076938 -3.9120230 7.677864 2.143157 -13.03161 -11.155251
## 8  4.513055 2.047693  2.0725428 6.703188 3.933333 -15.80328 -10.892349
## 9  5.178407 1.983756 -3.9120230 7.677864 2.208274 -15.80328 -11.381897
## 10 5.139322 2.063058  3.6095656 7.377759 4.163337 -12.89523 -10.785377
## 11 4.967728 2.008214 -3.9120230 7.200425 2.990034 -13.88058 -10.892349
## 12 3.916015 2.060514 -3.9120230 6.131226 2.472806 -14.93331 -11.512925
## 13 4.808111 2.010895  4.7621739 6.866933 4.672829 -12.93174 -10.892349
## 14 4.695925 2.018895  3.4011974 6.429719 3.934002 -14.14959 -11.258283
## 15 5.060694 2.041220  5.3303004 5.247024 4.704110 -13.40340 -10.346655
## 16 5.270946 2.013569 -0.1392621 7.146772 3.107850 -13.24453  -9.903488
## 17 5.141079 2.032088  4.5752263 6.927558 4.064862 -12.90725 -10.542147
## 18 5.015291 2.039921  3.6975914 7.098376 4.098021 -13.70218 -10.819778
## 19 4.581902 2.073172  4.9126549 6.476972 4.478278 -13.22772 -11.205441
## 20 4.374498 2.034706  4.1415462 6.476972 4.322603 -14.87106 -11.061850
## 21 4.774913 2.080691  3.1884166 6.253829 4.028095 -13.11738 -11.061850
## 22 4.632785 2.043814  4.5212450 6.291569 4.479401 -13.07357 -12.764689
## 23 4.745888 2.060514  3.1822118 6.672033 4.055124 -13.10256 -10.752120
## 24 4.827513 2.032088  6.9077553 6.040255 3.898134 -12.93588  -9.875872
## 25 4.305416 2.155245  4.1407509 6.109248 4.365588 -12.99573 -10.855405
## 26 4.097672 2.042518  3.7232809 5.991465 4.392179 -13.05470 -11.322305
## 27 4.332048 2.020222  3.3214324 6.214608 3.647207 -14.45607 -11.107460
## 28 4.198705 2.028148  4.7492705 6.086775 4.188276 -13.05470 -10.654264
## 29 4.239887 2.014903  4.1141472 6.131226 3.943330 -13.00015 -11.445267
## 30 4.461300 1.987874  3.3690185 6.345636 4.820532 -14.32634 -10.785377
## 31 4.919251 2.064328 -3.9120230 7.138867 2.580217 -14.12068 -12.764689
## 32 4.427239 2.041220 -3.9120230 6.606650 2.782952 -14.21004 -10.976432
## 33 4.112512 2.042518 -3.9120230 6.866933 1.798029 -13.86891 -10.976432
## 34 4.750136 2.096790  2.0175661 6.917706 3.767660 -12.93174 -11.258283
## 35 4.443827 2.081938  1.1052568 6.429719 3.673766 -12.84473 -10.976432
## 36 4.517431 2.064328 -3.9120230 6.672033 2.953935 -14.39891 -10.892349
## 37 4.091006 2.043814 -3.9120230 6.253829 3.506558 -15.11014 -11.754997
## 38 4.169761 2.063058 -3.9120230 6.345636 2.907440 -13.05939 -11.018229
## 39 4.320151 2.036012 -3.9120230 6.429719 2.471208 -13.23890 -11.849798
## 40 4.215086 2.051556 -3.9120230 6.363028 2.849778 -12.94422 -10.892349
## 41 3.949319 2.143589 -3.9120230 6.152733 2.625358 -13.16319 -14.152383
## 42 3.931826 2.104134 -3.9120230 6.194405 2.330756 -15.80328 -11.205441
## 43 4.418841 2.078191 -3.9120230 6.620073 3.224267 -12.94841 -12.206073
## 44 3.919991 2.028148 -3.9120230 6.152733 2.768515 -13.08314 -11.445267
## 45 5.196838 2.075684 -3.9120230 7.265430 2.644681 -14.03865 -10.625034
## 46 4.910447 2.038620  4.7405748 6.856462 4.499810 -13.47192 -10.654264
## 47 4.680278 2.048982  2.5610958 6.756932 3.642359 -13.66709 -10.688750
## 48 4.300003 2.013569 -3.9120230 6.646391 2.870044 -13.81551 -11.849798
## 49 4.565389 2.039921  4.6293749 6.756932 4.124881 -14.45607 -10.936312
## 50 4.527209 2.041220  3.6041382 6.363028 3.964469 -13.78595 -11.155251
## 51 4.647271 2.064328  1.4906544 6.507278 3.210844 -13.14768 -11.107460
## 52 4.237001 2.091864  1.8381658 6.309918 3.332205 -14.36196 -14.152383
## 53 4.289089 2.012233  4.2668963 6.214608 3.036554 -13.06882 -10.855405
## 54 4.262680 1.991976  3.4177267 6.152733 3.463541 -13.11738 -12.764689
##          NO2        NO3        Br         Cl         Fl        SO4        Na
## 1  -13.93092  -9.817310 -12.91942  -8.235041 -10.015537  -9.363492 -6.979470
## 2  -14.31309 -11.760106 -13.09280  -8.344998 -10.149388  -9.885648 -7.078470
## 3  -13.69329  -9.959000 -13.48621  -8.417488 -10.458613 -10.243165 -7.382206
## 4  -13.04078  -8.989071 -13.60040  -8.440737 -10.404363  -9.973910 -7.241213
## 7  -13.83776  -9.856604 -14.90022  -9.743071 -11.205441 -11.198115 -8.360685
## 8  -15.08136 -12.534577 -15.48682 -12.116232 -11.493123 -12.562748 -8.835946
## 9  -14.17935  -7.819967 -13.65000  -9.052858 -10.504968 -10.710924 -7.216567
## 10 -13.22772 -10.747458 -13.93317  -9.000971  -9.664471 -10.756803 -6.930678
## 11 -14.93637 -10.268771 -14.41153  -9.585761 -10.658510 -11.212821 -8.068873
## 12 -14.50866 -10.988197 -15.48682 -10.015537 -10.953310 -11.381897 -8.319154
## 13 -13.00903  -8.886570 -14.11932  -9.243897 -10.191170 -10.876349 -7.939656
## 14 -13.69329 -10.315977 -15.79601 -10.268771 -11.273909 -12.639937 -9.511445
## 15 -12.90725  -8.056369 -13.66709  -8.791788 -10.183201 -10.710924 -7.496422
## 16 -11.45466  -8.514038 -13.58440  -8.667017 -10.035877 -10.277454 -7.376421
## 17 -11.94988  -7.789553 -13.55315  -8.507312 -10.223693 -10.306955 -7.103282
## 18 -14.75456  -8.972789 -13.77629  -8.917395  -9.833961 -10.608707 -6.761265
## 19 -13.12236  -8.561866 -13.93317  -8.917186 -10.538366 -10.785377 -7.908529
## 20 -14.38644  -8.900667 -14.03740  -9.217365 -10.804890 -10.919599 -7.757987
## 21 -13.06882  -7.875745 -13.72933  -8.757489 -10.424364 -10.970601 -7.420468
## 22 -14.81247  -8.104240 -13.90544  -8.823467 -10.576832 -10.897740 -7.684699
## 23 -12.93588  -8.428715 -13.80556  -8.748917 -10.508624 -10.824791 -7.667536
## 24 -11.73357  -8.278432 -13.42347  -8.390292 -10.172675 -10.650036 -7.204263
## 25 -13.00015  -8.675515 -14.17076  -8.934377 -10.427736 -11.049191 -7.918280
## 26 -13.21119  -9.123337 -14.20707  -9.147300 -10.697561 -11.408565 -8.190274
## 27 -13.13241  -8.746267 -14.18948  -9.421061 -10.850237 -10.860600 -7.776376
## 28 -13.21119  -8.933383 -14.24476  -9.052858 -10.620927 -11.258283 -8.086106
## 29 -13.25020  -8.917738 -14.22649  -9.115239 -10.645825 -11.155251 -8.057136
## 30 -12.46744  -8.574171 -13.91976  -9.107829 -10.693146 -10.988197 -7.688587
## 31 -15.85573 -10.417652 -14.45796  -9.303553 -10.549751 -11.000102 -7.930663
## 32 -14.27755 -10.319003 -14.64105  -9.451139 -10.738198 -10.804890 -7.850380
## 33 -13.42347  -8.594868 -15.15111  -8.376106 -10.231992 -10.850237 -7.661326
## 34 -13.77629 -10.596635 -14.53295  -9.064977 -10.504968 -10.724468 -7.819272
## 35 -14.46752  -8.523810 -14.13572  -9.182054 -10.444772 -11.169336 -7.821006
## 36 -14.81247 -10.693146 -14.90022  -9.672376 -10.545942 -11.012150 -7.574122
## 37 -14.75456 -10.309953 -14.97706  -9.772459 -11.055501 -11.068240 -8.167939
## 38 -14.55366 -10.592643 -14.55785  -9.361163 -10.780558 -11.134489 -8.384587
## 39 -15.85573 -10.534599 -14.58554  -9.593066 -10.947612 -11.205441 -8.477747
## 40 -14.93637 -10.604667 -14.69982  -9.342730 -10.780558 -11.127663 -8.383826
## 41 -13.47192 -10.633299 -15.15111  -9.664471 -10.809828 -11.250561 -8.369463
## 42 -14.75456 -10.795086 -15.48682 -10.099502 -11.155251 -10.860600 -8.109847
## 43 -14.65022 -10.637457 -14.48299  -9.198372 -10.629158 -11.055501 -8.135051
## 44 -13.65851 -10.804890 -15.48682  -9.588677 -10.976432 -11.107460 -8.370965
## 45 -15.16258 -11.963911 -15.55848 -10.497695 -10.988197 -12.649240 -7.763909
## 46 -12.41433  -8.320989 -14.03740  -9.217365  -9.997798 -11.000102 -7.744840
## 47 -13.16841  -9.582854 -14.45796  -9.179045 -10.104380 -10.845096 -7.864545
## 48 -13.56865 -10.075463 -14.61180  -9.317926 -10.742817 -10.994132 -8.351429
## 49 -13.67575  -9.648845 -14.18948  -9.348454 -10.569020 -10.834892 -7.567454
## 50 -13.61666  -8.289480 -13.87845  -9.400291 -10.534599 -10.886987 -7.656624
## 51 -14.27755  -8.300913 -13.86575  -8.940603 -10.391248 -10.625034 -7.519088
## 52 -13.74785  -9.777736 -14.90022  -9.729534 -10.804890 -11.024345 -7.768115
## 53 -14.69982  -8.989977 -14.32467  -9.171418 -10.671358 -11.227747 -8.046346
## 54 -13.30269  -8.843859 -14.18948  -9.161169 -10.701995 -11.155251 -8.055215
##             K         Ca         Mg        TDN
## 1   -7.766051  -8.022238  -8.787134  -9.456241
## 2   -8.035378  -8.037332  -8.829236 -10.349775
## 3   -7.891573  -7.896082  -8.755051  -9.482149
## 4   -8.113042  -8.110599  -8.928128  -8.818298
## 7   -9.054748  -9.217365 -10.089817  -9.600424
## 8   -9.348454  -9.288302 -10.136681 -10.701995
## 9   -8.771460  -9.145096  -9.938079  -7.789645
## 10  -7.996265  -8.749182  -9.514152 -10.031321
## 11  -9.045818  -8.922421  -9.828380  -9.833961
## 12  -9.387078  -9.686765 -10.280365 -10.504968
## 13  -9.086107  -9.401501 -10.061312  -8.746606
## 14 -11.963911 -10.662775 -11.917891  -9.963238
## 15  -8.472746  -9.585761 -10.104380  -7.951879
## 16  -8.471646  -8.706074  -9.571310  -8.250990
## 17  -8.301493  -8.613313  -9.351904  -7.712952
## 18  -8.898417  -9.145564  -9.841452  -8.825078
## 19  -8.588147  -8.941879  -9.685156  -8.482792
## 20  -9.073313  -9.331379 -10.054310  -8.785073
## 21  -8.607972  -9.059570  -9.746484  -7.829059
## 22  -8.717060  -9.076448  -9.811820  -8.095199
## 23  -8.449292  -8.802498  -9.545813  -8.326573
## 24  -8.326965  -9.201351  -9.706277  -8.069307
## 25  -8.946546  -9.473005 -10.035877  -8.558015
## 26  -9.114512  -9.636519 -10.207299  -9.003326
## 27  -8.958205  -9.444798  -9.986869  -8.645027
## 28  -8.954838  -9.630412 -10.207299  -8.759265
## 29  -8.921904  -9.547213 -10.089817  -8.831904
## 30  -8.852666  -9.322390 -10.004413  -8.454218
## 31  -8.497768  -9.963238 -10.881654 -10.322038
## 32  -8.896355  -8.982966  -9.774215  -9.889585
## 33  -8.780566  -9.866192 -10.448215  -8.497391
## 34  -8.427544  -8.636236  -9.452412 -10.154516
## 35  -8.832308  -9.254292  -9.984698  -8.440232
## 36  -8.667847  -9.083672  -9.854697 -10.085009
## 37  -9.226470  -9.334770 -10.104380 -10.089817
## 38  -8.934050  -9.242864  -9.997798 -10.077841
## 39  -8.815283  -8.916455  -9.636519 -10.292096
## 40  -8.841816  -9.161836  -9.850895 -10.038162
## 41  -9.262687  -9.658191 -10.292096 -10.549751
## 42  -9.575624  -9.761988 -10.340443 -10.274551
## 43  -8.687970  -8.943601  -9.701363 -10.434516
## 44  -9.258481  -9.656627 -10.237563 -10.343544
## 45  -8.333543  -8.629517  -9.376395 -10.384754
## 46  -8.981740  -9.350753 -10.038162  -8.213392
## 47  -8.860388  -9.185160 -10.097072  -9.276480
## 48  -8.773105  -9.000209  -9.726179  -9.893537
## 49  -8.880598  -8.988542  -9.711216  -9.390664
## 50  -8.829858  -9.192726  -9.991226  -8.228262
## 51  -8.589931  -9.069275  -9.819146  -8.239561
## 52  -9.198847  -9.578510 -10.170061  -9.746484
## 53  -8.908009  -9.510095 -10.063656  -8.845697
## 54  -8.846253  -9.446063 -10.038162  -8.811564
biplot(pca2,scaling=1)

biplot(pca2,scaling=2)

# note different scaling factors, but solution remains same


##############################
# some follow-up suggestions #
# test PCA-axes for effect of landuse or stream size (as log(Q)) using ANOVA or ANCOVA
# correlate PCA-axes with other potential "controlling" variables (e.g. TDN, canopy cover) to give "meta-dimensions" more meaning
# useful function envfit() to relate additional variables to the ordination space

Redundancy analysis (RDA) in R

Two involved matrices: one dependent, one independent.

(Note: For correlation of two matrices see CCorA (Canonical correlation analysis) \(\neq\) CCA.)

Redundancy: The proportion of total variance of in the response variables that can be explained linear combinations of predictors.

Two steps:

  1. MLRs relate each response variable to the independent matrix and predict the response.
  2. The matrix of predicted response variables (same size as original: n objects * p variables) is subject to PCA.

The response variables are constrained to be linear combinations of the predictors first! The PCA can only “ordinate” variation of the responses that is relatable to predictors.

Significance tests for the overall model, for the various RDA-axes and for the individual predictors are available (permutation-based).

Site scores:

Redundancy analysis (RDA) in R

zwc<-scale(wc) # must at least be centered even if dimensionally homogeneous!
xmat<-data.frame(logQ=log(mara$Q),logTDN=log(mara$TDN),canopy=mara$canopy)[-c(5,6),]

rda<-rda(zwc~logQ+logTDN+canopy,data=xmat) # take care: confusing X and Y argument names

# actual RDA output check
summary(rda)
## 
## Call:
## rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          18.000     1.0000
## Constrained     6.217     0.3454
## Unconstrained  11.783     0.6546
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                         RDA1    RDA2    RDA3    PC1     PC2     PC3     PC4
## Eigenvalue            5.0745 0.94540 0.19681 4.9144 1.75906 1.00857 0.84980
## Proportion Explained  0.2819 0.05252 0.01093 0.2730 0.09773 0.05603 0.04721
## Cumulative Proportion 0.2819 0.33444 0.34537 0.6184 0.71612 0.77216 0.81937
##                           PC5     PC6     PC7     PC8     PC9     PC10     PC11
## Eigenvalue            0.70352 0.69192 0.52705 0.43646 0.29351 0.178515 0.115065
## Proportion Explained  0.03908 0.03844 0.02928 0.02425 0.01631 0.009917 0.006392
## Cumulative Proportion 0.85845 0.89689 0.92617 0.95042 0.96672 0.976642 0.983035
##                           PC12     PC13     PC14     PC15      PC16      PC17
## Eigenvalue            0.096526 0.080421 0.066249 0.041010 0.0137298 0.0074400
## Proportion Explained  0.005363 0.004468 0.003681 0.002278 0.0007628 0.0004133
## Cumulative Proportion 0.988397 0.992865 0.996546 0.998824 0.9995867 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1   RDA2    RDA3
## Eigenvalue            5.0745 0.9454 0.19681
## Proportion Explained  0.8163 0.1521 0.03166
## Cumulative Proportion 0.8163 0.9683 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.504409 
## 
## 
## Species scores
## 
##          RDA1      RDA2      RDA3        PC1        PC2        PC3
## cond -0.57654  0.612080 -0.021819 -6.413e-01 -2.623e-01 -2.670e-01
## pH    0.45481 -0.063413 -0.134308 -1.950e-01  5.990e-01 -4.274e-02
## turb -0.83125  0.001573  0.188277  3.949e-01  5.850e-01 -4.849e-01
## Alk   0.03483  0.465971 -0.022713 -7.910e-01 -5.520e-01 -3.544e-02
## TSS  -0.70702  0.219654  0.087775  4.023e-01  6.509e-01 -5.401e-01
## PO4  -0.31467 -0.045169 -0.308413  6.067e-02  9.303e-01  2.458e-01
## NH4  -0.53863  0.291317 -0.065288 -2.359e-01 -2.399e-01 -6.823e-01
## NO2  -0.80583 -0.193545  0.279175 -4.976e-02  3.768e-01 -1.452e-01
## NO3  -1.15497 -0.469987 -0.153789  4.834e-02  7.215e-02  2.538e-01
## Br   -0.92421  0.148249 -0.054180 -7.170e-01  1.920e-01 -2.714e-02
## Cl   -0.79416  0.002213  0.010841 -7.045e-01  3.496e-01  4.859e-01
## Fl   -0.72110  0.145664 -0.007607 -7.091e-01  4.037e-01  4.306e-02
## SO4  -0.49164  0.014945  0.139145 -9.938e-01  1.553e-01  3.314e-01
## Na   -0.70235  0.088094  0.050389 -9.311e-01  5.674e-02 -8.230e-02
## K    -0.32852  0.130777 -0.117953 -1.089e+00  1.871e-01 -7.620e-02
## Ca   -0.19361  0.501576 -0.081789 -1.047e+00 -6.644e-02 -2.463e-01
## Mg   -0.23802  0.339320 -0.048151 -1.079e+00  6.218e-02 -2.512e-01
## TDN  -1.22312 -0.408704 -0.142092 -4.770e-18  7.749e-17 -5.248e-17
## 
## 
## Site scores (weighted sums of species scores)
## 
##        RDA1     RDA2      RDA3       PC1      PC2        PC3
## 1  -0.89854  3.30672 -0.386506 -1.989464 -0.77705  0.7371045
## 2  -0.56847  4.04439  0.534394 -1.743560  0.92041 -0.8501656
## 3  -0.35983  3.03931 -0.034654 -1.545136 -1.28647  0.0427813
## 4  -0.77641  2.20741 -1.483632 -1.559658 -0.36335 -0.0291313
## 7   0.88348 -0.35314 -1.973257 -0.122169 -0.59904  0.1641181
## 8   1.89546  0.74944  2.486365  0.944714 -1.12313 -3.3055183
## 9  -0.52003  0.37308  0.198164 -0.477278 -2.80260  0.6379405
## 10 -0.63914  2.46873  0.151441 -0.963658  1.36394 -0.8199541
## 11  0.46308  1.27858 -1.064719 -0.391269 -0.86740 -0.4922014
## 12  1.44680 -1.16870  1.489863  0.041605 -0.33977 -0.3592235
## 13 -0.69258  0.08168  0.925485  0.819920  0.30530  0.3286962
## 14  1.41198 -2.32777  3.957086  3.261487 -0.87884  0.5640433
## 15 -1.23166 -0.83951  0.486820  0.463530  0.68386 -0.3343295
## 16 -1.26703  1.08821  0.262397 -0.531979 -0.42709  0.1677873
## 17 -1.58220  0.73379 -0.224828 -0.774707  0.18625 -0.6653627
## 18 -0.81300  1.11939 -0.295492 -0.233557  0.01539  0.1495812
## 19 -0.66646 -0.18332 -0.290191  0.319756  0.60538  0.1107652
## 20 -0.20474 -0.47474  1.516958  0.555638 -0.47108 -0.1279505
## 21 -0.93565 -0.55098 -1.434015 -0.025914  0.41162 -0.2975757
## 22 -0.56437 -0.70216 -1.466515  0.227191  0.38038  0.5203794
## 23 -0.83147  0.22256 -0.971790 -0.039558  0.27401 -0.0556500
## 24 -1.48587 -0.22424  0.743449 -0.289558  0.79299 -0.7436195
## 25 -0.37834 -1.18949 -1.029574  0.298653  1.42737 -0.4827891
## 26 -0.12196 -1.44806  0.564615  0.712294  0.85053 -0.1981804
## 27 -0.25191 -1.07885  1.937866  0.640716 -0.48858  0.2038573
## 28 -0.34372 -1.24137  0.490924  0.638891  0.55613 -0.3185752
## 29 -0.27388 -1.25221  0.563976  0.447446  0.47782  0.0534905
## 30 -0.72988 -0.59341  2.637512  0.500231 -0.27123 -0.4969936
## 31  0.90276 -0.02438 -0.951023 -0.098103 -0.24096  1.4932089
## 32  0.53785  0.31010 -0.000380 -0.195213 -0.44407  0.2430924
## 33  0.02874 -1.71644 -0.584641 -0.191690 -0.82265  1.3993386
## 34  0.10795  1.36339 -0.679776 -0.140368  1.14453  0.2044740
## 35 -0.24197 -0.66292 -2.716238  0.094133  0.30400 -0.1579230
## 36  0.67101  0.64767 -0.370543 -0.390286 -0.26097 -0.3706464
## 37  1.00078 -0.66312  1.406747  0.710440 -0.51072  0.5775779
## 38  0.77136 -0.28513 -1.624952  0.022634  0.37933  0.2749370
## 39  1.03318  0.20503 -2.567599 -0.054682 -0.06822  0.5670135
## 40  0.77241 -0.06114 -2.117887 -0.065082  0.24383  0.1888934
## 41  1.38747 -1.92709 -0.032796  0.002718  1.44558  1.2661054
## 42  1.49551 -1.30131  2.018641  0.082048 -0.65703 -0.4330951
## 43  0.71385  0.41829 -1.772975 -0.186986  0.86709  0.8222904
## 44  1.06214 -1.26145  0.642918  0.239372  0.27914  0.4771499
## 45  1.26992  2.30778 -2.592650  0.197937 -1.09176 -1.0417050
## 46 -0.94288 -0.05952  0.930304 -0.120353  0.25364 -1.2621277
## 47 -0.26421  0.23181  0.880294 -0.200536  0.30144 -0.4244519
## 48  0.52678 -0.11986  0.294289  0.069655 -0.25116  0.8955469
## 49 -0.22617  0.54673  1.879041 -0.101247  0.12653 -0.9326487
## 50 -0.55083 -0.68713 -0.007739  0.367544 -0.20260 -0.0379831
## 51 -0.54615 -0.33780 -2.156689 -0.035535 -0.10177  0.6969878
## 52  0.72450 -1.54996  2.095876  0.096846  0.56095  0.7707365
## 53 -0.05122 -0.96640 -1.060389  0.344503 -0.07383  0.0001981
## 54 -0.14642 -1.49251  0.796025  0.367647  0.26392  0.6797064
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         RDA1      RDA2      RDA3       PC1      PC2        PC3
## 1  -0.196294  0.827168  1.024793 -1.989464 -0.77705  0.7371045
## 2   0.456938  1.205241  1.334894 -1.743560  0.92041 -0.8501656
## 3  -0.155055  0.937232  0.556790 -1.545136 -1.28647  0.0427813
## 4  -0.350139 -0.151953  0.153012 -1.559658 -0.36335 -0.0291313
## 7   0.645108 -0.731268 -0.602499 -0.122169 -0.59904  0.1641181
## 8   1.399964 -0.094265 -0.277382  0.944714 -1.12313 -3.3055183
## 9  -0.949046 -0.756474 -1.035017 -0.477278 -2.80260  0.6379405
## 10  0.400897  1.170318 -0.584368 -0.963658  1.36394 -0.8199541
## 11  0.683071 -0.182817 -0.671894 -0.391269 -0.86740 -0.4922014
## 12  1.415316 -1.241293  1.341637  0.041605 -0.33977 -0.3592235
## 13 -0.900162  1.328320  0.066534  0.819920  0.30530  0.3286962
## 14  0.080004  1.535836  0.591574  3.261487 -0.87884  0.5640433
## 15 -1.179265 -0.030352  0.132810  0.463530  0.68386 -0.3343295
## 16 -1.033267  0.516938 -0.347991 -0.531979 -0.42709  0.1677873
## 17 -1.121367 -0.744455 -0.188840 -0.774707  0.18625 -0.6653627
## 18 -0.662565  0.751154  0.297991 -0.233557  0.01539  0.1495812
## 19 -0.920492  0.457760  0.672744  0.319756  0.60538  0.1107652
## 20 -0.586320  0.202346  0.905904  0.555638 -0.47108 -0.1279505
## 21 -0.911546 -0.594593 -1.513758 -0.025914  0.41162 -0.2975757
## 22 -0.758713 -0.515269 -0.966421  0.227191  0.38038  0.5203794
## 23 -0.928113  0.205313  0.237043 -0.039558  0.27401 -0.0556500
## 24 -0.903223 -0.450166 -0.113183 -0.289558  0.79299 -0.7436195
## 25 -0.368012 -0.496738 -0.746479  0.298653  1.42737 -0.4827891
## 26 -0.089848 -0.322454 -0.143368  0.712294  0.85053 -0.1981804
## 27 -0.697730  0.173662  0.847072  0.640716 -0.48858  0.2038573
## 28 -0.346198 -0.336436  0.185893  0.638891  0.55613 -0.3185752
## 29 -0.263150 -0.581729  0.747294  0.447446  0.47782  0.0534905
## 30 -0.791706 -0.061225  0.832429  0.500231 -0.27123 -0.4969936
## 31  0.923261  0.244046 -0.132397 -0.098103 -0.24096  1.4932089
## 32  0.493747  0.270364  0.091327 -0.195213 -0.44407  0.2430924
## 33 -0.228132 -1.320830 -0.056533 -0.191690 -0.82265  1.3993386
## 34  0.325305  1.783462 -0.729323 -0.140368  1.14453  0.2044740
## 35 -0.351369 -0.719932 -1.173188  0.094133  0.30400 -0.1579230
## 36  0.816178  0.094178 -0.641906 -0.390286 -0.26097 -0.3706464
## 37  0.562113  0.968540 -0.932109  0.710440 -0.51072  0.5775779
## 38  0.838366 -0.001746 -0.613313  0.022634  0.37933  0.2749370
## 39  1.028180  0.325781 -1.505883 -0.054682 -0.06822  0.5670135
## 40  0.826971 -0.078181 -0.607127 -0.065082  0.24383  0.1888934
## 41  1.462892 -0.923351  0.369467  0.002718  1.44558  1.2661054
## 42  1.320082 -1.062884 -0.162151  0.082048 -0.65703 -0.4330951
## 43  0.949646  0.830000 -1.161163 -0.186986  0.86709  0.8222904
## 44  1.047701 -0.638899  1.505043  0.239372  0.27914  0.4771499
## 45  0.483612  1.565549  0.343759  0.197937 -1.09176 -1.0417050
## 46 -0.529654 -1.122990 -0.261496 -0.120353  0.25364 -1.2621277
## 47 -0.008242 -0.330815  1.324407 -0.200536  0.30144 -0.4244519
## 48  0.420789  0.528903  0.006754  0.069655 -0.25116  0.8955469
## 49  0.094762  0.175222 -0.120623 -0.101247  0.12653 -0.9326487
## 50 -0.887653 -0.216708  0.330209  0.367544 -0.20260 -0.0379831
## 51 -0.864402 -0.088076 -0.151145 -0.035535 -0.10177  0.6969878
## 52  0.586747 -0.983296  1.772374  0.096846  0.56095  0.7707365
## 53 -0.136695 -0.620761 -0.120578  0.344503 -0.07383  0.0001981
## 54 -0.143291 -0.697375 -0.111622  0.367647  0.26392  0.6797064
## 
## 
## Biplot scores for constraining variables
## 
##           RDA1    RDA2    RDA3 PC1 PC2 PC3
## logQ    0.5159 -0.6932  0.5034   0   0   0
## logTDN -0.9427 -0.3150 -0.1095   0   0   0
## canopy  0.4907 -0.1427 -0.8596   0   0   0
?cca.object
RsquareAdj(rda) # redundancy statistic (fractional amount of variation of the response data matrix explained by constraints)
## $r.squared
## [1] 0.3453742
## 
## $adj.r.squared
## [1] 0.3044601

# hypothesis tests #
# testing the first axis (global test)
anova(rda)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance      F Pr(>F)    
## Model     3   6.2167 8.4414  0.001 ***
## Residual 48  11.7833                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda,first=TRUE)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance      F Pr(>F)    
## RDA1      1   5.0745 20.672  0.001 ***
## Residual 48  11.7833                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing all axes sequentially (preceding axes are taken as constraints)
anova(rda,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for rda under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## RDA1      1   5.0745 20.6715  0.001 ***
## RDA2      1   0.9454  3.8511  0.021 *  
## RDA3      1   0.1968  0.8017  0.546    
## Residual 48  11.7833                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing the individual terms=constraints
anova(rda,by="terms",model="direct",perm.max=9999,step=1000)  # tests terms sequentially, order matters!
## Permutation test for rda under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## logQ      1   1.8545  7.5544  0.002 ** 
## logTDN    1   3.7199 15.1532  0.001 ***
## canopy    1   0.6424  2.6168  0.035 *  
## Residual 48  11.7833                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda,by="margin",model="direct",perm.max=9999,step=1000) # tests each term in full model (like drop1() function)
## Permutation test for rda under direct model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## logQ      1   1.1757  4.7895  0.001 ***
## logTDN    1   2.7447 11.1808  0.001 ***
## canopy    1   0.6424  2.6168  0.030 *  
## Residual 48  11.7833                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

###################
# making triplots #

# again various types of scaling for the plotting step:
# scaling 1 "distance triplot"
# only angles between constraints and responses reflect their correlations (not angles among responses)
# distances among sites reflect their Euclidean distances

# scaling 2 "correlation triplot"
# all angles between constraints and responses reflect correlations 
# distances among sites do not reflect their Euclidean distances

# in both scaling types sites can be projected on constraints and on responses
# factor constraints are shown as centroids instead of arrows, projecting works identical 

# scaling 3 is compromise

# build an RDA scaling type 1 triplot
plot(rda,scaling=1)


(sites<-scores(rda,choices=c(1,2),display="sites",scaling=1)) 
##           RDA1         RDA2
## 1  -0.47708706  0.757822982
## 2  -0.30183628  0.926881198
## 3  -0.19105498  0.696540192
## 4  -0.41224315  0.505887724
## 7   0.46909247 -0.080931541
## 8   1.00641502  0.171753365
## 9  -0.27611373  0.085501680
## 10 -0.33935950  0.565775956
## 11  0.24587596  0.293021499
## 12  0.76819362 -0.267838709
## 13 -0.36773357  0.018719325
## 14  0.74970240 -0.533469771
## 15 -0.65396251 -0.192396774
## 16 -0.67274089  0.249392401
## 17 -0.84008293  0.168166947
## 18 -0.43167242  0.256537512
## 19 -0.35386382 -0.042012613
## 20 -0.10870767 -0.108800496
## 21 -0.49679038 -0.126272434
## 22 -0.29965548 -0.160919100
## 23 -0.44147805  0.051005608
## 24 -0.78893823 -0.051391214
## 25 -0.20088493 -0.272602952
## 26 -0.06475602 -0.331860547
## 27 -0.13375163 -0.247246487
## 28 -0.18250376 -0.284494012
## 29 -0.14541984 -0.286976426
## 30 -0.38753906 -0.135994789
## 31  0.47932829 -0.005588143
## 32  0.28557445  0.071067374
## 33  0.01525870 -0.393369015
## 34  0.05731908  0.312456891
## 35 -0.12847735 -0.151924742
## 36  0.35627838  0.148430402
## 37  0.53137430 -0.151970659
## 38  0.40955951 -0.065345506
## 39  0.54857945  0.046989184
## 40  0.41011812 -0.014011206
## 41  0.73669169 -0.441643953
## 42  0.79405601 -0.298230203
## 43  0.37902353  0.095862731
## 44  0.56395239 -0.289094360
## 45  0.67427788  0.528888976
## 46 -0.50063077 -0.013640487
## 47 -0.14028494  0.053125379
## 48  0.27970035 -0.027469107
## 49 -0.12008689  0.125297375
## 50 -0.29247055 -0.157474788
## 51 -0.28998363 -0.077414833
## 52  0.38467865 -0.355214393
## 53 -0.02719593 -0.221475935
## 54 -0.07774435 -0.342049506
## attr(,"const")
## [1] 5.504409

(lcs<-scores(rda,choices=c(1,2),display="lc",scaling=1)) # fitted/constrained site scores
##            RDA1         RDA2
## 1  -0.104224288  0.189567711
## 2   0.242615728  0.276213147
## 3  -0.082327802  0.214791744
## 4  -0.185909806 -0.034824199
## 7   0.342526429 -0.167589611
## 8   0.743324762 -0.021603400
## 9  -0.503905543 -0.173366219
## 10  0.212860509  0.268209766
## 11  0.362683168 -0.041897420
## 12  0.751476304 -0.284475586
## 13 -0.477950182  0.304419950
## 14  0.042478939  0.351977893
## 15 -0.626142621 -0.006955898
## 16 -0.548623320  0.118470198
## 17 -0.595401031 -0.170611798
## 18 -0.351795708  0.172146937
## 19 -0.488744614  0.104907983
## 20 -0.311312613  0.046373049
## 21 -0.483994623 -0.136266840
## 22 -0.402846352 -0.118087635
## 23 -0.492790671  0.047052871
## 24 -0.479575026 -0.103167648
## 25 -0.195399738 -0.113840754
## 26 -0.047705540 -0.073898965
## 27 -0.370466946  0.039799254
## 28 -0.183817338 -0.077103367
## 29 -0.139722215 -0.133318827
## 30 -0.420364221 -0.014031283
## 31  0.490214711  0.055929670
## 32  0.262159763  0.061961130
## 33 -0.121128764 -0.302703414
## 34  0.172723906  0.408727984
## 35 -0.186562979 -0.164991565
## 36  0.433357781  0.021583490
## 37  0.298459410  0.221966799
## 38  0.445138750 -0.000400045
## 39  0.545922368  0.074661431
## 40  0.439088339 -0.017917235
## 41  0.776737283 -0.211610578
## 42  0.700910966 -0.243588247
## 43  0.504223931  0.190216605
## 44  0.556287538 -0.146420841
## 45  0.256778620  0.358787439
## 46 -0.281224871 -0.257363258
## 47 -0.004376186 -0.075815134
## 48  0.223422164  0.121212240
## 49  0.050314794  0.040156741
## 50 -0.471308404 -0.049664526
## 51 -0.458962873 -0.020184905
## 52  0.311539183 -0.225348633
## 53 -0.072579521 -0.142264087
## 54 -0.076081553 -0.159822114
## attr(,"const")
## [1] 5.504409

(species<-scores(rda,choices=c(1,2),display="sp",scaling=1)*0.5)
##             RDA1         RDA2
## cond -0.54292631  1.335388807
## pH    0.42828677 -0.138350484
## turb -0.78277590  0.003432039
## Alk   0.03279901  1.016618867
## TSS  -0.66578993  0.479223111
## PO4  -0.29631964 -0.098546059
## NH4  -0.50722452  0.635572107
## NO2  -0.75884500 -0.422261198
## NO3  -1.08762377 -1.025380733
## Br   -0.87032206  0.323438441
## Cl   -0.74785116  0.004828477
## Fl   -0.67905469  0.317798293
## SO4  -0.46297404  0.032605655
## Na   -0.66140095  0.192196915
## K    -0.30936605  0.285318316
## Ca   -0.18232269  1.094299339
## Mg   -0.22414502  0.740302771
## TDN  -1.15180016 -0.891678553
## attr(,"const")
## [1] 5.504409

(constraints<-scores(rda,choices=c(1,2),display="bp",scaling=1)*2)
##              RDA1        RDA2
## logQ    0.5477956 -0.31771737
## logTDN -1.0011204 -0.14438929
## canopy  0.5210455 -0.06541799
## attr(,"const")
## [1] 5.504409

plot(sites,asp=1,pch=21,bg=landuse,ylim=c(-1.5,1.5))
Arrows(x0=0,y0=0,x1=constraints[,1],y1=constraints[,2],lwd=1.5,col="blue")
text(constraints[,1:2]*1.1,label=rownames(constraints),pos=4,cex=0.8,col="blue")

Arrows(x0=0,y0=0,x1=species[,1],y1=species[,2],lwd=1,arr.length=0)
text(species[,1:2]*1.1,label=rownames(species),pos=4,cex=0.6)

Correspondence Analysis (CA)

Based on one matrix (usually a species or community matrix).

Considers unimodal responses to (unknown) environmental variables.

An indirect GA, resulting gradients are synthetic environmental gradients.

The basis for CA is weighted averaging from environmental and species tables. If env exists, then this can be done to extract bioindicatory information:

\[ u^*=\frac{y_1x_1+y_2x_2+...+y_nx_n}{y_1+y_2+...+y_n} \] A species optimum \(u^*\) is computed as an abundance-weighted means of a specific environmental variable over all sites at which a specific species is present.

This approach works best when:

CA uses a two-way weighted averaging with a theoretical environmental variable iteratively in several steps:

  1. Take arbitrary site scores.
  2. Derive species scores by weighted average of sites scores, for species k (of m): \[ u_k=\sum_{i=1}^n{y_{ki}x_i}/\sum_{i=1}^n{y_{ki}} \]
  3. From the species scores new site scores can be derived, for site i (of n): \[ x_i=\sum_{k=1}^m{y_{ki}u_k}/\sum_{k=1}^n{y_{ki}} \]
  4. Rescaling (standardization) of site and species scores.
  5. Repeat 2-3 several times until stabilisation of site and species scores = first CA axis.
  6. Similar procedure to construct second CA axis (uncorrelated to first).

Canonical Correspondence Analysis (CCA)

Two involved matrices: one dependent, one independent.

In the reciprocal averaging of CA a constraint is included:

The result are axes which inform about species-site relationships, but which also have maximized correlation with linear combinations of (environmental) predictors.

Site scores:

data(varespec) # a R dataset on vegetation
data(varechem) # soil chemistry
head(varespec)
##    Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv Descflex Betupube
## 18     0.55    11.13     0.00     0.00    17.80     0.07     0.00        0
## 15     0.67     0.17     0.00     0.35    12.13     0.12     0.00        0
## 24     0.10     1.55     0.00     0.00    13.47     0.25     0.00        0
## 27     0.00    15.13     2.42     5.92    15.97     0.00     3.70        0
## 23     0.00    12.68     0.00     0.00    23.73     0.03     0.00        0
## 19     0.00     8.92     0.00     2.42    10.28     0.12     0.02        0
##    Vacculig Diphcomp Dicrsp Dicrfusc Dicrpoly Hylosple Pleuschr Polypili
## 18     1.60     2.07   0.00     1.62     0.00      0.0     4.67     0.02
## 15     0.00     0.00   0.33    10.92     0.02      0.0    37.75     0.02
## 24     0.00     0.00  23.43     0.00     1.68      0.0    32.92     0.00
## 27     1.12     0.00   0.00     3.63     0.00      6.7    58.07     0.00
## 23     0.00     0.00   0.00     3.42     0.02      0.0    19.42     0.02
## 19     0.00     0.00   0.00     0.32     0.02      0.0    21.03     0.02
##    Polyjuni Polycomm Pohlnuta Ptilcili Barbhatc Cladarbu Cladrang Cladstel
## 18     0.13     0.00     0.13     0.12     0.00    21.73    21.47     3.50
## 15     0.23     0.00     0.03     0.02     0.00    12.05     8.13     0.18
## 24     0.23     0.00     0.32     0.03     0.00     3.58     5.52     0.07
## 27     0.00     0.13     0.02     0.08     0.08     1.42     7.63     2.55
## 23     2.12     0.00     0.17     1.80     0.02     9.08     9.22     0.05
## 19     1.58     0.18     0.07     0.27     0.02     7.23     4.95    22.08
##    Cladunci Cladcocc Cladcorn Cladgrac Cladfimb Cladcris Cladchlo Cladbotr
## 18     0.30     0.18     0.23     0.25     0.25     0.23     0.00     0.00
## 15     2.65     0.13     0.18     0.23     0.25     1.23     0.00     0.00
## 24     8.93     0.00     0.20     0.48     0.00     0.07     0.10     0.02
## 27     0.15     0.00     0.38     0.12     0.10     0.03     0.00     0.02
## 23     0.73     0.08     1.42     0.50     0.17     1.78     0.05     0.05
## 19     0.25     0.10     0.25     0.18     0.10     0.12     0.05     0.02
##    Cladamau Cladsp Cetreric Cetrisla Flavniva Nepharct Stersp Peltapht Icmaeric
## 18     0.08   0.02     0.02     0.00     0.12     0.02   0.62     0.02        0
## 15     0.00   0.00     0.15     0.03     0.00     0.00   0.85     0.00        0
## 24     0.00   0.00     0.78     0.12     0.00     0.00   0.03     0.00        0
## 27     0.00   0.02     0.00     0.00     0.00     0.00   0.00     0.07        0
## 23     0.00   0.00     0.00     0.00     0.02     0.00   1.58     0.33        0
## 19     0.00   0.00     0.00     0.00     0.02     0.00   0.28     0.00        0
##    Cladcerv Claddefo Cladphyl
## 18        0     0.25        0
## 15        0     1.00        0
## 24        0     0.33        0
## 27        0     0.15        0
## 23        0     1.97        0
## 19        0     0.37        0
head(varechem)
##       N    P     K    Ca    Mg    S    Al   Fe    Mn   Zn  Mo Baresoil Humdepth
## 18 19.8 42.1 139.9 519.4  90.0 32.3  39.0 40.9  58.1  4.5 0.3     43.9      2.2
## 15 13.4 39.1 167.3 356.7  70.7 35.2  88.1 39.0  52.4  5.4 0.3     23.6      2.2
## 24 20.2 67.7 207.1 973.3 209.1 58.1 138.0 35.4  32.1 16.8 0.8     21.2      2.0
## 27 20.6 60.8 233.7 834.0 127.2 40.7  15.4  4.4 132.0 10.7 0.2     18.7      2.9
## 23 23.8 54.5 180.6 777.0 125.8 39.5  24.2  3.0  50.1  6.6 0.3     46.0      3.0
## 19 22.8 40.9 171.4 691.8 151.4 40.8 104.8 17.6  43.6  9.1 0.4     40.5      3.8
##     pH
## 18 2.7
## 15 2.8
## 24 3.0
## 27 2.8
## 23 2.7
## 19 2.7
apply(varespec,1,sum) # approximate 100 (total cover), "absolute" abundance data
##     18     15     24     27     23     19     22     16     28     13     14 
##  89.20  89.82  94.21 125.61  90.46  81.27 109.76  88.52 110.70 101.89  81.65 
##     20     25      7      5      6      3      4      2      9     12     10 
##  64.11  94.06 103.38  94.77 110.90 106.67  84.83 119.13 122.60 119.80 122.37 
##     11     21 
## 112.84  99.17

# correspondence analysis #
# run a CA just based on the species data (unconstrained!)
vare.ca<-cca(X=varespec) # function also used for CCA, but here only one matrix X is supplied

summary(vare.ca,scaling=1)
## 
## Call:
## cca(X = varespec) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total           2.083          1
## Unconstrained   2.083          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.5249 0.3568 0.2344 0.19546 0.17762 0.12156 0.11549
## Proportion Explained  0.2520 0.1713 0.1125 0.09383 0.08526 0.05835 0.05544
## Cumulative Proportion 0.2520 0.4233 0.5358 0.62962 0.71489 0.77324 0.82868
##                           CA8     CA9    CA10    CA11    CA12    CA13     CA14
## Eigenvalue            0.08894 0.07318 0.05752 0.04434 0.02546 0.01710 0.014896
## Proportion Explained  0.04269 0.03513 0.02761 0.02129 0.01222 0.00821 0.007151
## Cumulative Proportion 0.87137 0.90650 0.93411 0.95539 0.96762 0.97583 0.982978
##                           CA15     CA16     CA17     CA18     CA19      CA20
## Eigenvalue            0.010160 0.007830 0.006032 0.004008 0.002865 0.0019275
## Proportion Explained  0.004877 0.003759 0.002896 0.001924 0.001375 0.0009253
## Cumulative Proportion 0.987855 0.991614 0.994510 0.996434 0.997809 0.9987341
##                            CA21      CA22      CA23
## Eigenvalue            0.0018074 0.0005864 0.0002434
## Proportion Explained  0.0008676 0.0002815 0.0001168
## Cumulative Proportion 0.9996017 0.9998832 1.0000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##                 CA1       CA2      CA3       CA4        CA5       CA6
## Callvulg  0.0303167 -1.597460  0.11455 -2.894569  0.1376073  2.291129
## Empenigr  0.0751030  0.379305  0.39303  0.023675  0.8568729 -0.400964
## Rhodtome  1.1052309  1.499299  3.04284  0.120106  3.2324306 -0.283510
## Vaccmyrt  1.4614812  1.622935  2.72375  0.231688  0.4604556  0.712538
## Vaccviti  0.1468014  0.313436  0.14696  0.243505  0.6868371 -0.147815
## Pinusylv -0.4820096  0.588517 -0.36020 -0.127094  0.4064754  0.386604
## Descflex  1.5348239  1.218806  1.87562 -0.001340 -1.3136979 -0.070731
## Betupube  0.6694503  1.951826  3.84017  1.389423  7.5959115 -0.244478
## Vacculig -0.0830789 -1.629259  1.05063  0.802648 -0.3058811 -1.625341
## Diphcomp -0.5446464 -1.037570  0.52282  0.940275  0.3682126 -1.082929
## Dicrsp    1.8120408  0.360290 -4.92082  3.088562  1.3867372  0.157815
## Dicrfusc  1.2704743 -0.562978 -0.39718 -2.929542  0.3848272 -2.408710
## Dicrpoly  0.7248118  1.409347  0.80341  1.915549  4.5674148  1.295447
## Hylosple  2.0062408  1.743883  2.27549  0.928884 -3.7648428  2.254851
## Pleuschr  1.3102086  0.583036 -0.01004  0.137298 -1.1216144  0.200422
## Polypili -0.3805097 -1.243904  0.54593  1.477188 -0.7276341 -0.387641
## Polyjuni  1.0133795  0.099043 -2.24697  1.510641  0.7729714 -3.062378
## Polycomm  0.8468241  1.321773  1.13585  1.140723  2.6836594 -0.605038
## Pohlnuta -0.0136453  0.589290 -0.35542  0.135481  0.9369707  0.397246
## Ptilcili  0.4223631  1.598584  3.43474  1.400065  6.3209491  0.198935
## Barbhatc  0.5018348  2.119334  4.57303  1.693188  8.1101807  0.645995
## Cladarbu -0.1531729 -1.483884  0.20024  0.193680  0.0734141  0.358926
## Cladrang -0.5502561 -1.084008  0.40552  0.724060 -0.3357992 -0.335924
## Cladstel -1.4373146  1.077753 -0.44397 -0.375926 -0.2421525  0.004212
## Cladunci  0.8151727 -1.006186 -1.82587 -1.389523  1.6046713  3.675908
## Cladcocc -0.2133215 -0.584429 -0.21434 -0.567886 -0.0003788 -0.145303
## Cladcorn  0.2631227 -0.177858 -0.44464  0.272422  0.3992282 -0.306738
## Cladgrac  0.1956947 -0.311167 -0.23894  0.379013  0.4933026  0.037581
## Cladfimb  0.0009213 -0.161418  0.18463 -0.435908  0.4831233 -0.143751
## Cladcris  0.3373031 -0.470369 -0.05093 -0.823855  0.7182250  0.636140
## Cladchlo -0.6200021  1.207278  0.21889  0.426447  1.9506082  0.120722
## Cladbotr  0.5647242  1.047333  2.65330  0.907734  4.4946805  1.201655
## Cladamau -0.6598144 -1.512880  0.83251  1.577699 -0.0407227 -1.419139
## Cladsp   -0.8209003  0.476164 -0.49752 -0.998241 -0.2393208  0.390785
## Cetreric  0.2458192 -0.689228 -1.68427 -0.131681  0.7439412  2.374535
## Cetrisla -0.3465221  1.362693  0.85897  0.396752  2.7526968  0.396591
## Flavniva -1.4391907 -0.833589 -0.12919  0.007071 -1.4841375  2.956977
## Nepharct  1.6813309  0.199484 -4.33509  2.229917  0.9561223 -5.472858
## Stersp   -0.5172793 -2.280900  0.99775  2.377013 -0.8892757 -1.441228
## Peltapht  0.4035858 -0.043265  0.04538  0.711040  0.1824679 -0.841227
## Icmaeric  0.0378754 -2.419595  0.72135  0.361302 -0.3736424 -2.092136
## Cladcerv -0.9232858 -0.005233 -1.22058  0.305290 -0.8142627  0.414135
## Claddefo  0.5190399 -0.496632 -0.15271 -0.695927  0.9042143  0.909191
## Cladphyl -1.2836161  1.155872 -0.79912 -0.741170 -0.1608002  0.490526
## 
## 
## Site scores (weighted averages of species scores)
## 
##          CA1      CA2       CA3      CA4        CA5      CA6
## 18 -0.108122 -0.53705  0.229574  0.24412  0.1405624 -0.14253
## 15  0.697118 -0.14441 -0.031788 -0.21743 -0.2738522 -0.08146
## 24  0.987603  0.15042 -1.348447  0.80472  0.3095168  0.46773
## 27  0.851765  0.49901  0.443559  0.12277 -0.4814871  0.07589
## 23  0.359881 -0.05608  0.145813  0.15087  0.2405263 -0.17770
## 19  0.003545  0.37017  0.027760  0.06168 -0.1158930 -0.03413
## 22  0.860732 -0.11504  0.110869 -1.02169  0.0772348 -0.60530
## 16  0.636936 -0.33250  0.001120 -0.79797  0.0130769 -0.54049
## 28  1.279352  0.81557  0.670053  0.23137 -0.8929976  0.41783
## 13 -0.195009 -0.80564  0.117686 -0.58286 -0.0007212  0.53071
## 14  0.528532 -0.70420 -0.517771 -0.86836  0.5713441  0.91671
## 20  0.382866 -0.18686 -0.004789  0.10156  0.0458125  0.21087
## 25  0.990715  0.11967 -1.110040  0.44929  0.1885902 -0.70694
## 7  -0.264704 -1.06013  0.334900  0.45973 -0.0326631 -0.19945
## 5  -0.428410 -1.20765  0.374344  0.74970 -0.2596294 -0.30467
## 6  -0.330534 -0.77498  0.130760  0.22391  0.0632686  0.09060
## 3  -0.899601  0.12075 -0.075742  0.03842 -0.1489585 -0.12031
## 4  -0.770294 -0.35351 -0.033779 -0.01795 -0.3007839  0.44303
## 2  -0.992193  0.50319 -0.157505 -0.07070 -0.1065172 -0.09928
## 9  -0.937173  0.78688 -0.258119 -0.19377 -0.0343535 -0.01259
## 12 -0.726413  0.49163 -0.157235 -0.08698 -0.0105774 -0.02801
## 10 -1.002083  0.71239 -0.236526 -0.18643 -0.0231666 -0.04928
## 11 -0.322647 -0.03871 -0.001297  0.09029 -0.1481448  0.06934
## 21  0.259527  0.80746  1.124258  0.36083  1.5437866  0.07051
# summary(vare.ca,scaling=2)
# again two different types of scaling are possible for biplots

# scaling 1 (distances among sites matter)
# distances among sites approximate their chi^2 distance
# close sites have similar species abundances
# a site, which is near a specific species, has a high contribution of that species 

# scaling 2 (relationships among species matter)
# distances among species approximate their chi^2 distance
# close species have similar abundances across sites
# a species, which is near a specific site, is more likely to be found at that site

plot(vare.ca,scaling=1)

plot(vare.ca,scaling=2)

plot(vare.ca,scaling=3) # a compromise

# for any scaling take care when interpreting species close to origin:
# these are "everywhere" or have optimum right at the origin (i.e., optimum with regard to both axis shown)

# for more controlled plotting
species.scores<-scores(vare.ca,display="species",scaling=2)
site.scores<-scores(vare.ca,display="sites",scaling=2)

plot(site.scores,col="black",pch=21,xlim=c(-2,2),ylim=c(-2,2))
text(species.scores,col="red",label=names(varespec),cex=0.7)

# post-hoc fitting of an environmental variable
names(varechem)
##  [1] "N"        "P"        "K"        "Ca"       "Mg"       "S"       
##  [7] "Al"       "Fe"       "Mn"       "Zn"       "Mo"       "Baresoil"
## [13] "Humdepth" "pH"
(ef<-envfit(vare.ca,varechem[,12:13],permutations=1999))
## 
## ***VECTORS
## 
##               CA1      CA2     r2 Pr(>r)   
## Baresoil  0.97947 -0.20161 0.2533 0.0505 . 
## Humdepth  0.91602  0.40112 0.4524 0.0020 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1999
plot(ef)


###########################
# canonical correspondence analysis #
vare.cca<-cca(Y=varespec,X=varechem) # note strange terminology of X and Y in vegan (don´t ask)
vare.cca<-cca(varespec~.,varechem) # hypothesis tests need formula interface (don´t ask)

summary(vare.cca,scaling=1)
## 
## Call:
## cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.0832      1.000
## Constrained    1.4415      0.692
## Unconstrained  0.6417      0.308
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.16285 0.14213 0.11795 0.08903 0.07029
## Proportion Explained  0.2107 0.1401 0.07817 0.06823 0.05662 0.04274 0.03374
## Cumulative Proportion 0.2107 0.3507 0.42890 0.49713 0.55375 0.59649 0.63023
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.02801 0.01495 0.006382 0.004015 0.003139 0.002955
## Cumulative Proportion 0.65825 0.67319 0.679576 0.683592 0.686730 0.689685
##                          CCA14     CA1     CA2     CA3     CA4     CA5     CA6
## Eigenvalue            0.004733 0.19776 0.14193 0.10117 0.07079 0.05330 0.03330
## Proportion Explained  0.002272 0.09493 0.06813 0.04857 0.03398 0.02559 0.01598
## Cumulative Proportion 0.691958 0.78689 0.85502 0.90359 0.93757 0.96315 0.97914
##                            CA7      CA8      CA9
## Eigenvalue            0.018868 0.015104 0.009488
## Proportion Explained  0.009057 0.007251 0.004554
## Cumulative Proportion 0.988195 0.995446 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.1628 0.1421 0.11795 0.08903 0.07029
## Proportion Explained  0.3045 0.2024 0.1130 0.0986 0.08183 0.06176 0.04877
## Cumulative Proportion 0.3045 0.5069 0.6198 0.7184 0.80027 0.86203 0.91080
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.04049 0.02160 0.009223 0.005803 0.004536 0.004271
## Cumulative Proportion 0.95128 0.97288 0.982107 0.987910 0.992446 0.996716
##                          CCA14
## Eigenvalue            0.004733
## Proportion Explained  0.003284
## Cumulative Proportion 1.000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##              CCA1     CCA2     CCA3      CCA4     CCA5     CCA6
## Callvulg  0.11374 -1.73246  4.15753  1.844837  3.13741 -1.15626
## Empenigr -0.27373  0.14089  0.09036 -1.134550 -0.40226  0.03525
## Rhodtome -1.59033 -0.11156  0.19187 -2.490433 -0.62292 -1.73616
## Vaccmyrt -1.92827  0.56943  0.75259 -0.244263 -1.65624 -2.05452
## Vaccviti -0.23029  0.22315 -0.13141 -0.960947  0.24441  0.02995
## Pinusylv  0.36674  0.48934  0.55326 -0.726273  0.85051 -0.21227
## Descflex -2.17952  0.50020 -0.40165  1.608949 -1.38617  1.28224
## Betupube -1.07326 -0.41989 -0.20569 -6.388345 -0.62954 -5.60316
## Vacculig  0.77560 -2.19992 -0.93608  0.469588 -2.78966  1.04277
## Diphcomp  0.14991 -1.65300 -1.03898 -1.412059 -0.78833  2.08551
## Dicrsp   -1.28302  0.42862 -4.34136  0.691801  4.43282  1.30777
## Dicrfusc -0.75393 -0.76901  2.04376 -0.684761  0.32655  2.14058
## Dicrpoly -0.79564  0.14903 -2.01239 -3.186680  2.23820 -3.43647
## Hylosple -2.75940  1.46965  0.12345  3.602353 -2.66866 -0.74851
## Pleuschr -1.39625  0.62360 -0.02266  0.817213 -0.19077  0.06281
## Polypili  0.21763 -0.84392 -1.27708 -0.747467 -0.15333  0.16978
## Polyjuni -0.91607  0.38917 -0.87255 -0.891253 -1.78446  1.17847
## Polycomm -1.34974  0.59359 -0.58214 -2.854381 -1.19037 -2.60320
## Pohlnuta -0.01435  0.46778 -0.34834 -0.931562  1.23465 -0.32446
## Ptilcili -0.86964 -0.22650 -0.14520 -5.594842 -0.48392 -5.05263
## Barbhatc -1.04773 -0.42524 -0.29330 -6.830157 -0.50320 -6.88497
## Cladarbu  0.31928 -1.31813 -0.06534  0.138503 -0.11811 -0.26229
## Cladrang  0.57516 -1.14185 -0.60438  0.280955 -0.47617  0.10938
## Cladstel  1.36834  1.29986  0.20555  0.179762 -0.04827  0.09185
## Cladunci -0.34820  0.11797 -0.03422 -1.037583  2.65119 -0.48962
## Cladcocc  0.33121 -0.25213  0.31806 -0.205437  0.09828  0.41903
## Cladcorn -0.34025  0.12974 -0.22432 -0.686052 -0.31883  0.57211
## Cladgrac -0.16429 -0.34432 -0.39566 -0.533216  0.70217 -0.07237
## Cladfimb  0.03022 -0.16993  0.47734 -0.696051 -0.10470 -0.11656
## Cladcris -0.20689  0.02979  1.04812 -1.124294  0.40186 -0.43505
## Cladchlo  0.66964  1.02386 -0.68975 -1.528620  0.49217 -0.75368
## Cladbotr -1.02718 -0.35199  0.48348 -3.528217  0.63524 -4.23041
## Cladamau -0.02415 -2.15363 -1.80591 -1.323302 -1.02050  2.39498
## Cladsp    1.03577  0.72454  0.76099  0.741439  1.75911  0.41843
## Cetreric  0.09754 -0.07199 -1.05941  0.315233  2.75328 -0.58261
## Cetrisla  0.24027  0.64936 -0.12182 -2.346145  0.48511 -2.31098
## Flavniva  1.31684 -1.19677 -1.15320  5.202080  1.07346 -7.81575
## Nepharct -1.15139  0.36798 -1.38414 -0.153781 -3.31081  2.49381
## Stersp    0.18370 -2.37391 -2.38790 -0.009846 -1.07525  1.39790
## Peltapht -0.60047  0.31181  0.12300 -0.899163 -0.76856  0.65021
## Icmaeric  0.26085 -2.83828 -1.06550 -0.409685 -1.20472  1.06913
## Cladcerv  1.06877 -0.10889 -0.78377  3.250753  0.01418 -3.50019
## Claddefo -0.45498 -0.03870  0.60324 -1.497542  0.85219 -0.63272
## Cladphyl  1.51291  2.08492  0.04117 -0.268421  0.27480  0.48796
## 
## 
## Site scores (weighted averages of species scores)
## 
##        CCA1     CCA2      CCA3     CCA4       CCA5     CCA6
## 18  0.11823 -0.57251 -0.164982 -0.22892 -0.1940177  0.07213
## 15 -0.64276 -0.10649  0.169910  0.11432  0.0521025  0.23988
## 24 -0.84786  0.25736 -1.189184  0.14813  1.3580797  0.22853
## 27 -0.99432  0.35227  0.034639  0.28730 -0.4232940 -0.02911
## 23 -0.39622 -0.09941 -0.054725 -0.43892 -0.1038877  0.02098
## 19 -0.07306  0.38585  0.006695 -0.02930 -0.1896126 -0.02464
## 22 -0.72351 -0.26482  0.855780 -0.16216  0.0893297  0.55882
## 16 -0.50071 -0.42518  0.666714 -0.05991  0.1632133  0.51821
## 28 -1.48534  0.62159  0.100450  0.70953 -0.6209907 -0.35786
## 13  0.26728 -0.79352  0.904036  0.45978  0.6372512 -0.27314
## 14 -0.30232 -0.37451  0.439688 -0.39404  0.9278442  0.04663
## 20 -0.36984 -0.13664 -0.135727 -0.13735  0.0942856  0.03259
## 25 -0.85602  0.13551 -0.587776 -0.01017  0.3304823  0.65496
## 7   0.36936 -1.08951 -0.372699  0.05638 -0.4616061  0.05740
## 5   0.44060 -1.21454 -0.658393  0.16630 -0.4226872  0.15976
## 6   0.39218 -0.69770 -0.189710 -0.03141 -0.0990151 -0.05450
## 3   0.88634  0.21282 -0.085773  0.09809 -0.2111364  0.08974
## 4   0.77344 -0.30247 -0.083929  0.80863  0.1228687 -0.94716
## 2   0.93351  0.60859  0.004559  0.01574 -0.1379707  0.08149
## 9   0.86982  0.91296  0.096369 -0.05063  0.0005494  0.01469
## 12  0.67010  0.58561  0.034417 -0.09231 -0.0424654  0.05488
## 10  0.93439  0.83587  0.093851 -0.06296 -0.0540425  0.05003
## 11  0.30814  0.02923 -0.059108  0.09765 -0.0281526 -0.01160
## 21 -0.47641  0.23201  0.003915 -1.44448 -0.2880137 -1.21174
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        CCA1     CCA2      CCA3     CCA4      CCA5     CCA6
## 18 -0.28028 -0.71553 -0.198602 -0.35622 -0.016645  0.28042
## 15 -0.12605  0.26839  0.183427 -0.19962 -0.026309 -0.23569
## 24 -0.57190  0.13619 -1.113919  0.21486  1.130853  0.07844
## 27 -1.12491  0.26289 -0.227400  0.40474 -0.210924  0.14882
## 23 -0.52704  0.05792  0.103917 -0.34088 -0.098759  0.13090
## 19 -0.44851  0.54086  0.013494 -0.37832 -0.048521 -0.27997
## 22 -0.54244 -0.36270  0.612068 -0.02209  0.194629  0.66118
## 16 -0.09855 -0.62779  0.413117 -0.16871 -0.053130 -0.07504
## 28 -1.37258  0.59298  0.200796  0.71143 -0.478758 -0.19023
## 13  0.10953 -0.73196  1.049988  0.47179  0.604495 -0.16296
## 14 -0.09320  0.10867  0.313802 -0.33147  0.232443 -0.11453
## 20 -0.45423  0.04379 -0.082409 -0.42046  0.381970 -0.22780
## 25 -0.59995  0.15944 -0.222688 -0.02782 -0.388701  0.24252
## 7   0.91721 -1.04185 -0.323016  0.13738 -0.567910 -0.03541
## 5   0.06432 -1.09164 -0.636766  0.01507 -0.151542  0.29546
## 6   0.27735 -0.30740 -0.130894  0.02489 -0.019959  0.10057
## 3   0.63365  0.06729 -0.206032  0.05714 -0.365798 -0.04820
## 4   0.56735 -0.42871 -0.189592  0.87651  0.160886 -0.84789
## 2   1.01789  0.50232  0.038999  0.09780 -0.003433  0.21274
## 9   1.01611  0.86649 -0.006136 -0.04395  0.239962  0.19821
## 12  0.29646  0.12958  0.378874 -0.10628  0.044242  0.11423
## 10  0.73605  0.86077 -0.016804  0.04149 -0.158371  0.07950
## 11  0.39119  0.19766 -0.018368 -0.05333 -0.024357 -0.11579
## 21 -0.45499 -0.12585 -0.070007 -1.04926 -0.070611 -0.65096
## 
## 
## Biplot scores for constraining variables
## 
##              CCA1     CCA2      CCA3      CCA4      CCA5       CCA6
## N        -0.14766 -0.28570  0.002715  0.066860 -0.086965  0.0304388
## P        -0.21110  0.31268 -0.065374  0.180759  0.063227 -0.0363527
## K        -0.24254  0.16634  0.145204  0.180743  0.111771 -0.0586722
## Ca       -0.29655  0.22782 -0.015240  0.037048  0.105769  0.0129929
## Mg       -0.28817  0.18393 -0.057371  0.040679  0.170979 -0.0017181
## S        -0.01594  0.22454  0.059879  0.167562  0.205056 -0.0496189
## Al        0.50996 -0.02564  0.015177  0.147400  0.055261 -0.1004202
## Fe        0.43000 -0.04760 -0.016976  0.099139 -0.023974 -0.0332227
## Mn       -0.47852  0.12132  0.045621  0.109904 -0.047628  0.0538484
## Zn       -0.23723  0.18092 -0.112151  0.130337  0.212656 -0.0003564
## Mo        0.13523 -0.05582 -0.063359  0.122240  0.177367 -0.0935487
## Baresoil -0.35558 -0.13762  0.055249 -0.196250  0.057225 -0.1051508
## Humdepth -0.46157  0.10891  0.109612 -0.051173 -0.001117 -0.0153218
## pH        0.32935  0.04056 -0.131693  0.007887 -0.049995 -0.0176316
summary(vare.cca,scaling=2)
## 
## Call:
## cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.0832      1.000
## Constrained    1.4415      0.692
## Unconstrained  0.6417      0.308
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.16285 0.14213 0.11795 0.08903 0.07029
## Proportion Explained  0.2107 0.1401 0.07817 0.06823 0.05662 0.04274 0.03374
## Cumulative Proportion 0.2107 0.3507 0.42890 0.49713 0.55375 0.59649 0.63023
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.02801 0.01495 0.006382 0.004015 0.003139 0.002955
## Cumulative Proportion 0.65825 0.67319 0.679576 0.683592 0.686730 0.689685
##                          CCA14     CA1     CA2     CA3     CA4     CA5     CA6
## Eigenvalue            0.004733 0.19776 0.14193 0.10117 0.07079 0.05330 0.03330
## Proportion Explained  0.002272 0.09493 0.06813 0.04857 0.03398 0.02559 0.01598
## Cumulative Proportion 0.691958 0.78689 0.85502 0.90359 0.93757 0.96315 0.97914
##                            CA7      CA8      CA9
## Eigenvalue            0.018868 0.015104 0.009488
## Proportion Explained  0.009057 0.007251 0.004554
## Cumulative Proportion 0.988195 0.995446 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.1628 0.1421 0.11795 0.08903 0.07029
## Proportion Explained  0.3045 0.2024 0.1130 0.0986 0.08183 0.06176 0.04877
## Cumulative Proportion 0.3045 0.5069 0.6198 0.7184 0.80027 0.86203 0.91080
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.04049 0.02160 0.009223 0.005803 0.004536 0.004271
## Cumulative Proportion 0.95128 0.97288 0.982107 0.987910 0.992446 0.996716
##                          CCA14
## Eigenvalue            0.004733
## Proportion Explained  0.003284
## Cumulative Proportion 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CCA1     CCA2      CCA3      CCA4      CCA5      CCA6
## Callvulg  0.075347 -0.93581  1.677742  0.695507  1.077518 -0.345001
## Empenigr -0.181340  0.07610  0.036462 -0.427727 -0.138153  0.010517
## Rhodtome -1.053549 -0.06026  0.077428 -0.938897 -0.213938 -0.518031
## Vaccmyrt -1.277428  0.30759  0.303704 -0.092088 -0.568820 -0.613023
## Vaccviti -0.152563  0.12054 -0.053031 -0.362279  0.083942  0.008938
## Pinusylv  0.242956  0.26432  0.223265 -0.273806  0.292102 -0.063335
## Descflex -1.443872  0.27019 -0.162082  0.606576 -0.476067  0.382590
## Betupube -0.711004 -0.22681 -0.083007 -2.408417 -0.216212 -1.671857
## Vacculig  0.513817 -1.18831 -0.377748  0.177035 -0.958084  0.311138
## Diphcomp  0.099310 -0.89289 -0.419273 -0.532348 -0.270745  0.622270
## Dicrsp   -0.849964  0.23153 -1.751924  0.260810  1.522412  0.390210
## Dicrfusc -0.499460 -0.41539  0.824743 -0.258156  0.112149  0.638702
## Dicrpoly -0.527090  0.08050 -0.812083 -1.201383  0.768689 -1.025365
## Hylosple -1.828026  0.79385  0.049816  1.358093 -0.916528 -0.223338
## Pleuschr -0.924978  0.33684 -0.009146  0.308091 -0.065518  0.018741
## Polypili  0.144172 -0.45586 -0.515356 -0.281796 -0.052660  0.050659
## Polyjuni -0.606869  0.21021 -0.352109 -0.336004 -0.612858  0.351629
## Polycomm -0.894165  0.32063 -0.234919 -1.076106 -0.408823 -0.776736
## Pohlnuta -0.009508  0.25268 -0.140571 -0.351201  0.424031 -0.096811
## Ptilcili -0.576115 -0.12234 -0.058593 -2.109265 -0.166198 -1.507591
## Barbhatc -0.694092 -0.22970 -0.118360 -2.574980 -0.172821 -2.054320
## Cladarbu  0.211517 -0.71201 -0.026366  0.052216 -0.040564 -0.078262
## Cladrang  0.381030 -0.61678 -0.243893  0.105921 -0.163536  0.032637
## Cladstel  0.906486  0.70213  0.082949  0.067771 -0.016579  0.027407
## Cladunci -0.230671  0.06372 -0.013810 -0.391170  0.910527 -0.146092
## Cladcocc  0.219419 -0.13619  0.128350 -0.077450  0.033754  0.125028
## Cladcorn -0.225404  0.07008 -0.090524 -0.258643 -0.109501  0.170706
## Cladgrac -0.108836 -0.18599 -0.159664 -0.201023  0.241156 -0.021594
## Cladfimb  0.020022 -0.09179  0.192626 -0.262413 -0.035959 -0.034780
## Cladcris -0.137056  0.01609  0.422960 -0.423861  0.138016 -0.129810
## Cladchlo  0.443621  0.55305 -0.278345 -0.576292  0.169030 -0.224882
## Cladbotr -0.680481 -0.19013  0.195105 -1.330144  0.218169 -1.262258
## Cladamau -0.015996 -1.16331 -0.728763 -0.498887 -0.350481  0.714608
## Cladsp    0.686166  0.39137  0.307091  0.279524  0.604150  0.124850
## Cetreric  0.064619 -0.03889 -0.427516  0.118844  0.945590 -0.173838
## Cetrisla  0.159171  0.35076 -0.049161 -0.884501  0.166607 -0.689545
## Flavniva  0.872373 -0.64645 -0.465365  1.961193  0.368671 -2.332045
## Nepharct -0.762768  0.19877 -0.558560 -0.057976 -1.137069  0.744096
## Stersp    0.121697 -1.28229 -0.963619 -0.003712 -0.369284  0.417103
## Peltapht -0.397796  0.16843  0.049634 -0.338986 -0.263955  0.194009
## Icmaeric  0.172805 -1.53313 -0.429975 -0.154452 -0.413750  0.319003
## Cladcerv  0.708032 -0.05882 -0.316283  1.225539  0.004871 -1.044377
## Claddefo -0.301412 -0.02090  0.243431 -0.564576  0.292677 -0.188788
## Cladphyl  1.002262  1.12620  0.016613 -0.101195  0.094379  0.145598
## 
## 
## Site scores (weighted averages of species scores)
## 
##       CCA1     CCA2      CCA3     CCA4     CCA5     CCA6
## 18  0.1785 -1.05988 -0.408835 -0.60721 -0.56492  0.24175
## 15 -0.9702 -0.19714  0.421046  0.30324  0.15171  0.80394
## 24 -1.2798  0.47645 -2.946863  0.39292  3.95433  0.76592
## 27 -1.5009  0.65216  0.085837  0.76207 -1.23251 -0.09756
## 23 -0.5981 -0.18404 -0.135611 -1.16425 -0.30249  0.07033
## 19 -0.1103  0.71431  0.016591 -0.07773 -0.55210 -0.08258
## 22 -1.0921 -0.49026  2.120668 -0.43014  0.26010  1.87287
## 16 -0.7558 -0.78712  1.652152 -0.15892  0.47523  1.73677
## 28 -2.2421  1.15075  0.248921  1.88204 -1.80814 -1.19935
## 13  0.4035 -1.46904  2.240249  1.21956  1.85549 -0.91541
## 14 -0.4563 -0.69333  1.089571 -1.04519  2.70161  0.15628
## 20 -0.5583 -0.25296 -0.336340 -0.36433  0.27453  0.10923
## 25 -1.2922  0.25087 -1.456542 -0.02698  0.96227  2.19508
## 7   0.5576 -2.01700 -0.923568  0.14954 -1.34406  0.19237
## 5   0.6651 -2.24847 -1.631533  0.44110 -1.23074  0.53544
## 6   0.5920 -1.29165 -0.470112 -0.08331 -0.28830 -0.18265
## 3   1.3379  0.39399 -0.212551  0.26020 -0.61477  0.30075
## 4   1.1675 -0.55997 -0.207980  2.14490  0.35776 -3.17436
## 2   1.4091  1.12669  0.011297  0.04175 -0.40173  0.27311
## 9   1.3130  1.69016  0.238808 -0.13429  0.00160  0.04923
## 12  1.0115  1.08413  0.085287 -0.24485 -0.12365  0.18392
## 10  1.4105  1.54744  0.232569 -0.16699 -0.15736  0.16768
## 11  0.4651  0.05411 -0.146473  0.25902 -0.08197 -0.03886
## 21 -0.7191  0.42952  0.009702 -3.83149 -0.83861 -4.06109
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        CCA1     CCA2     CCA3     CCA4      CCA5    CCA6
## 18 -0.42308 -1.32466 -0.49215 -0.94489 -0.048464  0.9398
## 15 -0.19026  0.49687  0.45454 -0.52951 -0.076603 -0.7899
## 24 -0.86328  0.25213 -2.76035  0.56993  3.292710  0.2629
## 27 -1.69805  0.48669 -0.56351  1.07358 -0.614147  0.4988
## 23 -0.79557  0.10723  0.25751 -0.90419 -0.287557  0.4387
## 19 -0.67702  1.00130  0.03344 -1.00351 -0.141279 -0.9383
## 22 -0.81881 -0.67147  1.51674 -0.05858  0.566703  2.2159
## 16 -0.14877 -1.16222  1.02373 -0.44751 -0.154699 -0.2515
## 28 -2.07190  1.09778  0.49758  1.88707 -1.394002 -0.6375
## 13  0.16534 -1.35508  2.60193  1.25142  1.760111 -0.5461
## 14 -0.14069  0.20118  0.77762 -0.87922  0.676806 -0.3838
## 20 -0.68566  0.08107 -0.20421 -1.11529  1.112185 -0.7635
## 25 -0.90562  0.29517 -0.55183 -0.07379 -1.131782  0.8128
## 7   1.38453 -1.92877 -0.80045  0.36440 -1.653585 -0.1187
## 5   0.09709 -2.02095 -1.57794  0.03999 -0.441247  0.9902
## 6   0.41866 -0.56908 -0.32436  0.06603 -0.058116  0.3371
## 3   0.95649  0.12458 -0.51056  0.15157 -1.065096 -0.1616
## 4   0.85641 -0.79366 -0.46982  2.32495  0.468453 -2.8417
## 2   1.53650  0.92994  0.09664  0.25941 -0.009995  0.7130
## 9   1.53381  1.60412 -0.01520 -0.11658  0.698700  0.6643
## 12  0.44751  0.23990  0.93887 -0.28191  0.128819  0.3828
## 10  1.11107  1.59354 -0.04164  0.11005 -0.461130  0.2664
## 11  0.59050  0.36592 -0.04552 -0.14145 -0.070919 -0.3881
## 21 -0.68681 -0.23299 -0.17348 -2.78317 -0.205599 -2.1817
## 
## 
## Biplot scores for constraining variables
## 
##              CCA1     CCA2      CCA3     CCA4      CCA5      CCA6
## N        -0.22290 -0.52891  0.006729  0.17735 -0.253216  0.102014
## P        -0.31866  0.57886 -0.162001  0.47947  0.184099 -0.121835
## K        -0.36612  0.30794  0.359824  0.47942  0.325444 -0.196637
## Ca       -0.44764  0.42176 -0.037765  0.09827  0.307969  0.043545
## Mg       -0.43499  0.34051 -0.142169  0.10790  0.497841 -0.005758
## S        -0.02406  0.41570  0.148384  0.44446  0.597063 -0.166296
## Al        0.76978 -0.04747  0.037610  0.39098  0.160905 -0.336554
## Fe        0.64909 -0.08811 -0.042067  0.26297 -0.069806 -0.111345
## Mn       -0.72232  0.22460  0.113052  0.29152 -0.138680  0.180471
## Zn       -0.35810  0.33493 -0.277916  0.34572  0.619191 -0.001195
## Mo        0.20413 -0.10334 -0.157007  0.32424  0.516439 -0.313525
## Baresoil -0.53675 -0.25477  0.136910 -0.52055  0.166621 -0.352409
## Humdepth -0.69673  0.20163  0.271625 -0.13574 -0.003252 -0.051350
## pH        0.49716  0.07509 -0.326341  0.02092 -0.145569 -0.059091
# again two different types of scaling are possible for triplots

# hypothesis tests #
# testing the first axis (global test)
anova(vare.cca)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## Model    14   1.44148 1.4441  0.034 *
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vare.cca,first=TRUE)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## CCA1      1   0.43887 6.1551  0.027 *
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing all axes sequentially (preceding axes are taken as constraints)
anova(vare.cca,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for cca under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## CCA1      1   0.43887 6.1551  0.035 *
## CCA2      1   0.29178 4.0921  0.405  
## CCA3      1   0.16285 2.2839  0.962  
## CCA4      1   0.14213 1.9934  0.983  
## CCA5      1   0.11795 1.6543  0.995  
## CCA6      1   0.08903 1.2486  1.000  
## CCA7      1   0.07029 0.9859  1.000  
## CCA8      1   0.05836 0.8185  1.000  
## CCA9      1   0.03114 0.4367  1.000  
## CCA10     1   0.01329 0.1865  1.000  
## CCA11     1   0.00836 0.1173  1.000  
## CCA12     1   0.00654 0.0917  1.000  
## CCA13     1   0.00616 0.0863  1.000  
## CCA14     1   0.00473 0.0664  1.000  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing the individual terms=constraints
anova(vare.cca,by="terms",model="direct",perm.max=9999,step=1000)  # tests terms sequentially, order matters!
## Permutation test for cca under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)   
## N         1   0.13001 1.8234  0.066 . 
## P         1   0.17996 2.5240  0.013 * 
## K         1   0.13643 1.9134  0.077 . 
## Ca        1   0.08614 1.2081  0.310   
## Mg        1   0.06760 0.9481  0.482   
## S         1   0.17545 2.4606  0.010 **
## Al        1   0.14797 2.0753  0.021 * 
## Fe        1   0.04982 0.6987  0.720   
## Mn        1   0.05080 0.7124  0.698   
## Zn        1   0.07167 1.0052  0.425   
## Mo        1   0.09585 1.3442  0.215   
## Baresoil  1   0.08683 1.2178  0.283   
## Humdepth  1   0.09502 1.3327  0.216   
## pH        1   0.06794 0.9529  0.485   
## Residual  9   0.64171                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vare.cca,by="margin",model="direct",perm.max=9999,step=1000) # tests each term in full model (like drop1() function)
## Permutation test for cca under direct model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## N         1   0.06302 0.8838  0.545  
## P         1   0.10327 1.4484  0.161  
## K         1   0.10385 1.4565  0.117  
## Ca        1   0.06169 0.8652  0.585  
## Mg        1   0.08830 1.2384  0.257  
## S         1   0.10388 1.4569  0.139  
## Al        1   0.06876 0.9644  0.455  
## Fe        1   0.04765 0.6683  0.777  
## Mn        1   0.08110 1.1375  0.308  
## Zn        1   0.06448 0.9043  0.550  
## Mo        1   0.07761 1.0885  0.361  
## Baresoil  1   0.08966 1.2574  0.280  
## Humdepth  1   0.13548 1.9001  0.041 *
## pH        1   0.06794 0.9529  0.476  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# quite a lot of variables in the constraining matrix, maybe selection would be adequate
# --> function ordistep()

###################
# making triplots #

# again various types of scaling for the plotting step:
# scaling 1 "distance triplot"
# sites can be projected on constraints
# sites close to centroid of factor constraint are more likely to possess the specific state (factor level)
# distances among sites reflect their Chi^2 distances

# scaling 2 "correlation triplot"
# species can be projected on constraints (to give their optimum)
# species close to centroid of factor constraint are more likely to be found in the respective sites
# distances among sites do not reflect their Chi^2 distances

# scaling 3 is compromise

plot(vare.cca,scaling=1,display=c("species","sites"))
plot(vare.cca,scaling=1,display=c("species","sites","bp"))

plot(vare.cca,scaling=2)

plot(vare.cca,scaling=3) # a compromise

# for any scaling take care when interpreting species close to origin:
# these are "everywhere" or have optimum right at the origin (i.e., optimum with regard to both axis shown)

# for more controlled plotting compute scores individually

(species<-scores(vare.cca,display="species",scaling=2))
##                  CCA1        CCA2
## Callvulg  0.075346828 -0.93580864
## Empenigr -0.181339634  0.07610159
## Rhodtome -1.053549303 -0.06026078
## Vaccmyrt -1.277428382  0.30758571
## Vaccviti -0.152563159  0.12053851
## Pinusylv  0.242955730  0.26432438
## Descflex -1.443872268  0.27018775
## Betupube -0.711004196 -0.22680891
## Vacculig  0.513817276 -1.18831302
## Diphcomp  0.099309627 -0.89289046
## Dicrsp   -0.849963650  0.23152562
## Dicrfusc -0.499459808 -0.41538963
## Dicrpoly -0.527089621  0.08050066
## Hylosple -1.828025570  0.79385044
## Pleuschr -0.924977860  0.33684271
## Polypili  0.144171903 -0.45585592
## Polyjuni -0.606869478  0.21021273
## Polycomm -0.894165298  0.32063347
## Pohlnuta -0.009508326  0.25267888
## Ptilcili -0.576115080 -0.12234488
## Barbhatc -0.694091803 -0.22969682
## Cladarbu  0.211516504 -0.71200604
## Cladrang  0.381029817 -0.61678185
## Cladstel  0.906485693  0.70213449
## Cladunci -0.230670950  0.06372124
## Cladcocc  0.219418695 -0.13619040
## Cladcorn -0.225403730  0.07008025
## Cladgrac -0.108836318 -0.18598627
## Cladfimb  0.020022402 -0.09178764
## Cladcris -0.137056217  0.01608948
## Cladchlo  0.443621232  0.55305019
## Cladbotr -0.680480562 -0.19013376
## Cladamau -0.015995645 -1.16331056
## Cladsp    0.686166357  0.39136838
## Cetreric  0.064619228 -0.03888872
## Cetrisla  0.159171117  0.35076227
## Flavniva  0.872372656 -0.64644751
## Nepharct -0.762767923  0.19876924
## Stersp    0.121696721 -1.28229477
## Peltapht -0.397796035  0.16842583
## Icmaeric  0.172805490 -1.53313439
## Cladcerv  0.708032101 -0.05881737
## Claddefo -0.301411810 -0.02090192
## Cladphyl  1.002261821  1.12619548

(lcs<-scores(vare.cca,display="lc",scaling=2)) # fitted site scores
##           CCA1        CCA2
## 18 -0.42308410 -1.32465581
## 15 -0.19026462  0.49687345
## 24 -0.86328447  0.25212611
## 27 -1.69805066  0.48669129
## 23 -0.79556975  0.10723446
## 19 -0.67702233  1.00129602
## 22 -0.81881442 -0.67146524
## 16 -0.14876518 -1.16222054
## 28 -2.07190438  1.09778134
## 13  0.16533766 -1.35507570
## 14 -0.14068713  0.20117513
## 20 -0.68566132  0.08107181
## 25 -0.90561895  0.29516717
## 7   1.38453057 -1.92876582
## 5   0.09708596 -2.02094655
## 6   0.41865544 -0.56908384
## 3   0.95649196  0.12457945
## 4   0.85641341 -0.79366326
## 2   1.53649728  0.92993610
## 9   1.53380971  1.60412304
## 12  0.44750585  0.23989837
## 10  1.11106621  1.59354430
## 11  0.59049933  0.36592318
## 21 -0.68681020 -0.23299449

(sites<-scores(vare.cca,display="sites",scaling=2)) # unfitted site scores
##          CCA1        CCA2
## 18  0.1784733 -1.05988423
## 15 -0.9702382 -0.19713866
## 24 -1.2798478  0.47644978
## 27 -1.5009195  0.65215591
## 23 -0.5980933 -0.18403623
## 19 -0.1102881  0.71431421
## 22 -1.0921288 -0.49026245
## 16 -0.7558244 -0.78712482
## 28 -2.2421137  1.15075172
## 13  0.4034539 -1.46904155
## 14 -0.4563466 -0.69332916
## 20 -0.5582740 -0.25295922
## 25 -1.2921584  0.25086578
## 7   0.5575516 -2.01700439
## 5   0.6650822 -2.24846553
## 6   0.5919915 -1.29165187
## 3   1.3379214  0.39399435
## 4   1.1674982 -0.55996655
## 2   1.4091299  1.12668910
## 9   1.3129824  1.69015916
## 12  1.0115068  1.08413112
## 10  1.4104517  1.54743675
## 11  0.4651334  0.05410696
## 21 -0.7191331  0.42951771

(constraints<-scores(vare.cca,choices=c(1,2),display="bp",scaling=2))
##                 CCA1        CCA2
## N        -0.22289549 -0.52891359
## P        -0.31865680  0.57886019
## K        -0.36611630  0.30794063
## Ca       -0.44763820  0.42176452
## Mg       -0.43499067  0.34051494
## S        -0.02405594  0.41569739
## Al        0.76977830 -0.04747443
## Fe        0.64908911 -0.08811470
## Mn       -0.72232495  0.22459941
## Zn       -0.35809796  0.33492792
## Mo        0.20413234 -0.10334347
## Baresoil -0.53674724 -0.25477013
## Humdepth -0.69673051  0.20163085
## pH        0.49715734  0.07508792

plot(sites,col="black",pch=21,xlim=c(-2,2),ylim=c(-2,2))
text(species,col="red",label=names(varespec),cex=0.7)
Arrows(x0=0,y0=0,x1=constraints[,1],y1=constraints[,2],lwd=1.5,col="blue")
text(constraints[,1:2]*1.1,label=rownames(constraints),pos=4,cex=0.8,col="blue")

Methods based on distance/dissimilarity/similarity

Why use dissimilarity-based methods?

Whenever compositional turnover at regional scale is of interest!

Diversity partitioning into \(\alpha\), \(\beta\), \(\gamma\):

\(\beta\)-diversity is linked to species-accumulation curves.

(Dis)similarity coefficients

Euclidean and Bray-Curtis (Steinhaus, Sorenson) are two examples of many dissimilarity indices, here computed between any two sites for \(p\) variables:

\[ d_{Euc}=\sqrt{\Delta{x_1}^2+\Delta{x_2}^2+...+\Delta{x_p}^2} \] where \(\Delta{x}\) is the distance between two sites along any \(X\)-variable.


\[ d_{BC}=\frac{\sum_{j=1}^p{|x_{1j}-x_{2j}|}}{\sum_{j=1}^p{(x_{1j}+x_{2j})}} \] which includes a notable normalization of differences between two sites with regard to any variable used.

(Dis)similarity coefficients: some properties


\[ D_{(a,b)}+D_{(b,c)}\geq{D_{(a,c)}} \] (think of a triangle: the distance from A to C via B must be larger than the direct distance from A to C)

\[ D_{(a,b)}+D_{(b,c)}<{D_{(a,c)}} \] … in which case perfect projections into Euclidean space are not possible. Two often used transformations to make semimetric coefficients analyzable in Euclidean space are adding a constant or \(\sqrt(D)\).

(Dis)similarity coefficients


Using dissimilarity: Cluster analysis

Hierarchical group-forming (clustering) based on pairwise (dis)similarity.

Three types:

  1. Divisive: start with all, successively split into 2.
  2. Agglomerative: start with individual observations and cluster pairwise, continue grouping clusters.
  3. Non-hierarchical: e.g. K-means clustering (forcing K clusters)


Various “linkage rules” to group clusters in agglomerative clustering:

lipids<-read.table(file="data/BacterialMembrane.txt",header=TRUE)
names(lipids)
##  [1] "case"        "isolate"     "temperature" "FA1_SAnb"    "FA2_MU"     
##  [6] "FA3_SAnb"    "FA4_SAb"     "FA5_SAb"     "FA6_SAnb"    "FA7_MU"     
## [11] "FA8_SAb"     "FA9_SAnb"    "FA10_MU"     "FA11_SAb"    "FA12_SAnb"  
## [16] "FA13_MUb"    "FA14_MU"     "FA15_SAnb"   "FA16_MU"     "FA17_MU"    
## [21] "sum_MU"      "sum_SA"      "sum_SAbran"  "sum_SAnbran" "SA_branprop"
lip_data<-lipids[,4:20]
# abbreviations: MU = mono-unsaturated, SA = saturated, nb = non-branched, b = branched
# SA_branprop = proportion of branched and saturated FAs
# theory: unsaturated and branched FAs increase fluidity of membrane
# test (i) adaptation vs. (ii) acclimatization to temperature by changing FA composition of membranes

isolate<-factor(lipids$isolate)
temperature<-factor(lipids$temperature)
combifac<-factor(paste(temperature,"_",isolate,sep=""))

## distance matrix: Euclidean distance on arcsine-sqr-data
as_lip_data<-asin(sqrt(lip_data)) # very old-school, better don´t do ;-)
lipids_distE<-vegdist(as_lip_data, method="euclidean") 

## alternatively: Bray-Curtis distance with proportional data
lipids_distBC<-vegdist(lip_data, method="bray")

# cluster analysis #
# various agglomeration methods available and the choice is important, explore method=
# "single": nearest neighbour counts (good for gradients, but makes chains)
# "complete": all group members must be close, farthest group member counts (makes small spheric groups, good to find outliers)
# "average": compromise average strategy, new member joins at mean distance to all group members, actually UPGMA
# "ward.D": aims at minimizing within-group sums of squares of distances
lipids_cluster<-hclust(lipids_distBC, method = "ward.D")
lipids_cluster$height<-sqrt(lipids_cluster$height) # may help

# then compare effect of agglomeration method on dendrogram
plot(lipids_cluster, hang = -1, labels = combifac, ylab = "BC")


# which agglomeration method (and thus dendrogram) is best?
cophenetic(lipids_cluster) # linkage distances in dendrogram
##             1          2          3          4          5          6          7
## 2  0.16296148                                                                  
## 3  0.40441248 0.40441248                                                       
## 4  0.40441248 0.40441248 0.31840951                                            
## 5  1.40874876 1.40874876 1.40874876 1.40874876                                 
## 6  1.40874876 1.40874876 1.40874876 1.40874876 0.20777591                      
## 7  1.40874876 1.40874876 1.40874876 1.40874876 0.27324836 0.27324836           
## 8  0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 9  0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 10 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 11 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 12 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 13 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 14 1.40874876 1.40874876 1.40874876 1.40874876 0.49516799 0.49516799 0.49516799
## 15 1.40874876 1.40874876 1.40874876 1.40874876 0.49516799 0.49516799 0.49516799
## 16 1.40874876 1.40874876 1.40874876 1.40874876 0.49516799 0.49516799 0.49516799
## 17 0.51155960 0.51155960 0.51155960 0.51155960 1.40874876 1.40874876 1.40874876
## 18 0.51155960 0.51155960 0.51155960 0.51155960 1.40874876 1.40874876 1.40874876
## 19 0.51155960 0.51155960 0.51155960 0.51155960 1.40874876 1.40874876 1.40874876
## 20 0.40441248 0.40441248 0.25690465 0.31840951 1.40874876 1.40874876 1.40874876
## 21 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834 0.56164834 0.56164834
## 22 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 23 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 24 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 25 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876 1.40874876 1.40874876
## 26 1.40874876 1.40874876 1.40874876 1.40874876 0.39214779 0.39214779 0.39214779
## 27 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834 0.56164834 0.56164834
## 28 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834 0.56164834 0.56164834
##             8          9         10         11         12         13         14
## 2                                                                              
## 3                                                                              
## 4                                                                              
## 5                                                                              
## 6                                                                              
## 7                                                                              
## 8                                                                              
## 9  0.23904776                                                                  
## 10 0.29578401 0.29578401                                                       
## 11 0.23904776 0.20171061 0.29578401                                            
## 12 0.05572996 0.23904776 0.29578401 0.23904776                                 
## 13 0.23904776 0.04176804 0.29578401 0.20171061 0.23904776                      
## 14 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876           
## 15 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 0.16915116
## 16 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 0.24387639
## 17 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 18 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 19 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 20 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 21 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834
## 22 0.49724975 0.49724975 0.49724975 0.49724975 0.49724975 0.49724975 1.40874876
## 23 0.49724975 0.49724975 0.49724975 0.49724975 0.49724975 0.49724975 1.40874876
## 24 0.49724975 0.49724975 0.49724975 0.49724975 0.49724975 0.49724975 1.40874876
## 25 0.29578401 0.29578401 0.18736413 0.29578401 0.29578401 0.29578401 1.40874876
## 26 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 0.49516799
## 27 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834
## 28 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834
##            15         16         17         18         19         20         21
## 2                                                                              
## 3                                                                              
## 4                                                                              
## 5                                                                              
## 6                                                                              
## 7                                                                              
## 8                                                                              
## 9                                                                              
## 10                                                                             
## 11                                                                             
## 12                                                                             
## 13                                                                             
## 14                                                                             
## 15                                                                             
## 16 0.24387639                                                                  
## 17 1.40874876 1.40874876                                                       
## 18 1.40874876 1.40874876 0.12548408                                            
## 19 1.40874876 1.40874876 0.30920822 0.30920822                                 
## 20 1.40874876 1.40874876 0.51155960 0.51155960 0.51155960                      
## 21 0.56164834 0.56164834 1.40874876 1.40874876 1.40874876 1.40874876           
## 22 1.40874876 1.40874876 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 23 1.40874876 1.40874876 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 24 1.40874876 1.40874876 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 25 1.40874876 1.40874876 0.97985304 0.97985304 0.97985304 0.97985304 1.40874876
## 26 0.49516799 0.49516799 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834
## 27 0.56164834 0.56164834 1.40874876 1.40874876 1.40874876 1.40874876 0.46298482
## 28 0.56164834 0.56164834 1.40874876 1.40874876 1.40874876 1.40874876 0.46298482
##            22         23         24         25         26         27
## 2                                                                   
## 3                                                                   
## 4                                                                   
## 5                                                                   
## 6                                                                   
## 7                                                                   
## 8                                                                   
## 9                                                                   
## 10                                                                  
## 11                                                                  
## 12                                                                  
## 13                                                                  
## 14                                                                  
## 15                                                                  
## 16                                                                  
## 17                                                                  
## 18                                                                  
## 19                                                                  
## 20                                                                  
## 21                                                                  
## 22                                                                  
## 23 0.06318240                                                       
## 24 0.25154043 0.25154043                                            
## 25 0.49724975 0.49724975 0.49724975                                 
## 26 1.40874876 1.40874876 1.40874876 1.40874876                      
## 27 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834           
## 28 1.40874876 1.40874876 1.40874876 1.40874876 0.56164834 0.23758016
plot(lipids_distBC,cophenetic(lipids_cluster))

cor(lipids_distBC,cophenetic(lipids_cluster),method="spearman")
## [1] 0.8600969

cutree(lipids_cluster,k=4)
##  [1] 1 1 1 1 2 2 2 3 3 3 3 3 3 2 2 2 1 1 1 1 4 3 3 3 3 2 4 4
cutree(lipids_cluster,h=0.05)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12  9 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27

plot(lipids_cluster, hang = -1, labels = combifac, ylab = "BC")
rect.hclust(lipids_cluster,k=4) # or specify height h instead of k

Using dissimilarity: Principal Coordinate Analysis (PCoA)

(aka metric scaling)

Technique to project sites onto Cartesian (Euclidean) coordinates from pairwise dissimilarities.

Imagine a road distance matrix between major towns. A 2D-projection of all pairwise distances produces a map (which gets better with less mountains and straighter roads in the area).

PCoA can be regarded as the inverse operation to computing a distance matrix. However, it can´t produce more “axes” than (n-1) as this is the maximum number of dimensions needed to completely represent n objects. A dissimilarity matrix and a consecutive PCoA are efficient means to condense data with p>>n.

PCoA is an important step in hypothesis-testing methods in the distance domain (see later).


PCoA is based on an eigenvector decomposition like PCA:

# using the same distance matrix as cluster analysis 
pcoa<-cmdscale(lipids_distBC,k=2,eig=TRUE,add=TRUE)
# argument add=TRUE means a constant is added to distances to avoid negative eigenvalues

cumsum(pcoa$eig/sum(pcoa$eig)) # contributions of various PCoA axes
##  [1] 0.5108797 0.6713554 0.7243661 0.7626502 0.7943804 0.8232193 0.8499307
##  [8] 0.8685434 0.8850589 0.8969014 0.9077931 0.9172892 0.9256833 0.9338634
## [15] 0.9419715 0.9496211 0.9562653 0.9622121 0.9677456 0.9731763 0.9782572
## [22] 0.9829655 0.9875741 0.9920859 0.9960663 1.0000000 1.0000000 1.0000000
# first two axes cover 67% of variation of distances

pcoa$points # the site scores (coordinates in reduced space)
##               [,1]         [,2]
##  [1,]  0.039034416 -0.063039145
##  [2,]  0.052535090 -0.034222468
##  [3,]  0.027399434 -0.087663904
##  [4,] -0.009878077 -0.004018133
##  [5,] -0.212039581  0.030589392
##  [6,] -0.216084698  0.020983290
##  [7,] -0.259894649  0.051819890
##  [8,]  0.165391223  0.116835950
##  [9,]  0.150166208  0.118147706
## [10,]  0.185980907  0.073277130
## [11,]  0.128185960  0.114502815
## [12,]  0.163993524  0.117098404
## [13,]  0.149963125  0.117191012
## [14,] -0.219431512  0.064195917
## [15,] -0.210959910  0.057130086
## [16,] -0.251459219  0.070724865
## [17,]  0.040185127 -0.182894184
## [18,]  0.034791299 -0.174443966
## [19,]  0.036014921 -0.211286639
## [20,]  0.012082138 -0.103177494
## [21,] -0.160218882 -0.026898255
## [22,]  0.202366734 -0.008237486
## [23,]  0.198633241 -0.005774575
## [24,]  0.198199603 -0.068765798
## [25,]  0.194183841  0.040168498
## [26,] -0.233079249  0.023071464
## [27,] -0.078445944 -0.018582519
## [28,] -0.127615070 -0.026731854

col.isolate<-isolate
levels(col.isolate)<-c("white","red")
col.isolate<-as.character(col.isolate)

pch.temperature<-as.numeric(as.character(temperature))
pch.temperature[pch.temperature==6]<-21
pch.temperature[pch.temperature==28]<-23

plot(pcoa$points,pch=pch.temperature,bg=col.isolate)
legend("topleft",pch=c(21,21,23,23),pt.bg=c("white","red","white","red"),
       legend=c("6°C - warm isolate","6°C - cold isolate","28°C - warm isolate","28°C - cold isolate"),cex=0.6)

# how to relate species (=fatty acids) to ordination?
wascores(pcoa$points,lip_data) # as weighted averages of site (=sample) scores
text(wascores(pcoa$points,lip_data),labels=names(lip_data),cex=0.5)

ordisurf(pcoa$points,lip_data$FA7_MU,col="darkgreen",add=TRUE) # as contourplot

plot(envfit(pcoa$points,lip_data)) # take care: behaviour of species not necessarily monotonous in ordination space

Using dissimilarity: Non-metric multidimensional scaling (NMDS)

An iterative search for an ordination configuration in a low-dimensional space with inter-object distances representing the observed dissimilarity/distance matrix as well as possible.

The measure of fit used is not actual dissimarity/distance but rank order of dissimilarities.

Iterative procedure:

  1. Compute matrix of dissimilarities.
  2. Decide on k, the number of dimensions.
  3. Arrange objects in a random starting configuration.
  4. Compute a measure of fit that expresses the match between inter-object distances of the configuration and the observed dissimilarities. A Shepard-plot shows residuals as stress which is inversely related to fit. The measure of fit is computed using ranks of observed dissimilarities and configuration distances.
  5. Reiteratively reposition the objects in the low-dimensional space and recompute fit to improve the match between inter-object distances and observed dissimilarities.
  6. A final configuration is achieved when no more repositioning improves the fit. Steps 3) to 6) may be repeated with different random starting positions to avoid getting trapped in local minima.
# non-metric multidimensional scaling (NMDS) #
# note that all resulting axes are equally important, plots may be rotated as needed

# using the same distance matrix as cluster analysis 
#mds_lipids<-isoMDS(lipids_distBC, k = 2, maxit=1000)
mds_lipids<-metaMDS(comm=lip_data, distance="bray", k = 2,trymax=100)
## Run 0 stress 0.06823601 
## Run 1 stress 0.06291168 
## ... New best solution
## ... Procrustes: rmse 0.05036073  max resid 0.1893593 
## Run 2 stress 0.07362809 
## Run 3 stress 0.06291168 
## ... Procrustes: rmse 7.14281e-06  max resid 2.47217e-05 
## ... Similar to previous best
## Run 4 stress 0.0682362 
## Run 5 stress 0.08036014 
## Run 6 stress 0.07362809 
## Run 7 stress 0.06291168 
## ... Procrustes: rmse 4.44794e-07  max resid 1.629583e-06 
## ... Similar to previous best
## Run 8 stress 0.07228103 
## Run 9 stress 0.06823697 
## Run 10 stress 0.06796483 
## Run 11 stress 0.06823625 
## Run 12 stress 0.07362809 
## Run 13 stress 0.06291168 
## ... Procrustes: rmse 1.326353e-05  max resid 3.172257e-05 
## ... Similar to previous best
## Run 14 stress 0.07228097 
## Run 15 stress 0.06291168 
## ... New best solution
## ... Procrustes: rmse 3.589165e-06  max resid 7.27529e-06 
## ... Similar to previous best
## Run 16 stress 0.06823585 
## Run 17 stress 0.06796483 
## Run 18 stress 0.06917428 
## Run 19 stress 0.06823585 
## Run 20 stress 0.06291168 
## ... Procrustes: rmse 4.989635e-06  max resid 1.47675e-05 
## ... Similar to previous best
## *** Solution reached

mds_lipids$stress # approximately 6.8% of dissimilarities remain unrepresented
## [1] 0.06291168
mds_lipids$points # the site scores
##           MDS1         MDS2
## 1  -0.07281623 -0.042093146
## 2  -0.08834359 -0.016851191
## 3  -0.04897736 -0.075177960
## 4   0.01313624  0.070365419
## 5   0.31339301  0.068577462
## 6   0.31363381  0.088897350
## 7   0.38016017  0.110505209
## 8  -0.24831670  0.131109175
## 9  -0.23159160  0.129025762
## 10 -0.27432580  0.085572655
## 11 -0.20184629  0.128763888
## 12 -0.24619739  0.129699656
## 13 -0.23159160  0.129025757
## 14  0.34897958 -0.018654271
## 15  0.33013239 -0.013738077
## 16  0.40718523  0.009552504
## 17 -0.09652212 -0.204858333
## 18 -0.09263457 -0.201461350
## 19 -0.10685476 -0.291406627
## 20 -0.03121302 -0.120771275
## 21  0.28955893 -0.193173036
## 22 -0.27468239 -0.003772281
## 23 -0.27256597 -0.003232314
## 24 -0.27317859 -0.065288401
## 25 -0.26414048  0.056903940
## 26  0.34983139  0.211431835
## 27  0.11512929 -0.043948014
## 28  0.19465842 -0.055004337
## attr(,"centre")
## [1] TRUE
## attr(,"pc")
## [1] TRUE
## attr(,"halfchange")
## [1] TRUE
## attr(,"internalscaling")
## [1] 3.663216
stressplot(mds_lipids,lipids_distBC)

(gof<-goodness(mds_lipids,statistic="distance")) # goodness of fit by NMDS for each sample
##  [1] 0.011062354 0.013697609 0.007607018 0.017971020 0.012576325 0.014203934
##  [7] 0.008087433 0.009154932 0.009137032 0.008889532 0.011453452 0.010106705
## [13] 0.009137032 0.010449041 0.011982104 0.013634174 0.011813176 0.011758968
## [19] 0.017537587 0.016442166 0.017081825 0.005768584 0.005300513 0.008717898
## [25] 0.007322299 0.013383859 0.011959534 0.013389617

# plotting of MDS scores
plot(mds_lipids$points,asp=1,xlab="MDS dimension 1",ylab="MDS dimension 2")
points(mds_lipids$points,pch=pch.temperature,bg=col.isolate,cex=gof*200) # to check goodness of fit, large symbol means bad fit

legend("topleft",pch=c(21,21,23,23),pt.bg=c("white","red","white","red"),
    legend=c("6°C - warm isolate","6°C - cold isolate","28°C - warm isolate","28°C - cold isolate"),cex=0.7)

# relating underlying variables to ordination in nMDS
wascores(mds_lipids$points,lip_data) # as weighted averages of site (=sample) scores
text(wascores(pcoa$points,lip_data),labels=names(lip_data),cex=0.5)

ordisurf(mds_lipids,lip_data$FA7_MU,col="darkgreen",add=TRUE)

plot(envfit(pcoa$points,lip_data)) # take care: behaviour of species not necessarily monotonous in ordination space

# some more useful graphical tools
plot(mds_lipids$points,pch=pch.temperature,bg=col.isolate,asp=1,xlab="MDS dimension 1",ylab="MDS dimension 2",cex=2)
ordispider(mds_lipids$points,groups=combifac)
ordihull(mds_lipids$points,groups=combifac)
ordiellipse(mds_lipids$points,groups=combifac)
ordicluster(mds_lipids$points,cluster=lipids_cluster)
plot(mds_lipids$points,pch=pch.temperature,bg=col.isolate,asp=1,xlab="MDS dimension 1",ylab="MDS dimension 2",cex=2)
ordispider(mds_lipids$points,groups=combifac)

Hypothesis tests in the distance world

All known study designs with factors or continuous predictors may be transferred to the distance domain.

Permutational MANOVA (PERMANOVA)

Uses a test statistic based on distances within groups (to a group centroid or averaged among all pairs) versus distances from group centroids to the overall centroid.


  1. The within-group sum of squares is the sum of squared distances from individual replicates to their group centroid.
  2. The among-group sum of squares is the sum of squared distances from group centroids to the overall centroid.
  3. A (pseudo-) F-value is computed using number of groups (a) and the total number of observations (N) as: \[ F=\frac{SS_A/(a-1)}{SS_w/(N-a)} \]
  4. Significance is assessed by recomputing the test-statistic after permutations of group assignment.


Prerequisite similar to ANOVA: homogeneous dispersion (multivariate variance, cloud shape). Tested using within-group distances to centroids.

Permutational MANOVA (PERMANOVA)

# PERMANOVA - non-parametric permutational MANOVA #
# a multivariate hypothesis test: 
# two factors and 1 multivariate response = "membrane FA composition"

# first testing for homogeneity of dispersion (homogeneous distances to group centroids)
disp.check<-betadisper(lipids_distBC,combifac)
disp.check$distances
##  [1] 0.13346669 0.11147756 0.15422004 0.15102228 0.06074253 0.06835717
##  [7] 0.05569824 0.03298576 0.01986574 0.05190567 0.02237701 0.03134896
## [13] 0.01877437 0.05215781 0.04500705 0.06406574 0.10040564 0.10078771
## [19] 0.13388533 0.10118002 0.14105434 0.06286176 0.06364213 0.03720424
## [25] 0.09455924 0.15755988 0.06427489 0.02935132
boxplot(disp.check$distances~combifac)

anova(lm(disp.check$distances~combifac))
## Analysis of Variance Table
## 
## Response: disp.check$distances
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## combifac   3 0.004967 0.0016556  0.8283 0.4913
## Residuals 24 0.047972 0.0019988
anova(disp.check)
## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## Groups     3 0.004967 0.0016556  0.8283 0.4913
## Residuals 24 0.047972 0.0019988
permutest(disp.check)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     3 0.004967 0.0016556 0.8283    999  0.505
## Residuals 24 0.047972 0.0019988

# the actual PERMANOVA
adonis(lipids_distBC~isolate*temperature)
## 
## Call:
## adonis(formula = lipids_distBC ~ isolate * temperature) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                     Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## isolate              1   0.06599 0.06599   7.988 0.09987  0.004 ** 
## temperature          1   0.38394 0.38394  46.476 0.58105  0.001 ***
## isolate:temperature  1   0.01258 0.01258   1.522 0.01903  0.193    
## Residuals           24   0.19826 0.00826         0.30005           
## Total               27   0.66077                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Canonical analysis of principal coordinates (CAP)

Following scheme taken from description of db-RDA. Essentially 3 steps:

  1. Computation of a (square) dissimilarity matrix D from (Cartesian) raw data. Choose coefficient well!
  2. PCoA based on D recreates Cartesian coordinates with a dimensionality imposed by n (or p if p<n). If D was semimetric, (minor) axes with negative eigenvalues may occur.
  3. RDA on the PCoA-axes and constraints of choice (factors/dummy variables, continuous predictors). Instead of RDA other constrained methods are possible as well, e.g. discriminant analysis, other classification routines, etc.

Canonical analysis of principal coordinates (CAP)

cap<-capscale(lipids_distBC~isolate*temperature)
summary(cap)
## 
## Call:
## capscale(formula = lipids_distBC ~ isolate * temperature) 
## 
## Partitioning of squared Bray distance:
##               Inertia Proportion
## Total          0.7018     1.0000
## Constrained    0.4670     0.6654
## Unconstrained  0.2348     0.3346
## 
## Eigenvalues, and their contribution to the squared Bray distance 
## 
## Importance of components:
##                         CAP1    CAP2     CAP3   MDS1    MDS2    MDS3    MDS4
## Eigenvalue            0.4027 0.05493 0.009286 0.1233 0.03707 0.02196 0.01708
## Proportion Explained  0.5739 0.07827 0.013231 0.1756 0.05283 0.03129 0.02434
## Cumulative Proportion 0.5739 0.65214 0.665369 0.8410 0.89384 0.92514 0.94947
##                          MDS5     MDS6     MDS7     MDS8     MDS9    MDS10
## Eigenvalue            0.01462 0.007081 0.004262 0.002948 0.001993 0.001796
## Proportion Explained  0.02083 0.010090 0.006073 0.004200 0.002840 0.002558
## Cumulative Proportion 0.97030 0.980392 0.986465 0.990665 0.993505 0.996063
##                          MDS11     MDS12     MDS13     MDS14     MDS15
## Eigenvalue            0.001280 0.0008798 0.0003408 0.0002224 3.548e-05
## Proportion Explained  0.001824 0.0012536 0.0004856 0.0003169 5.055e-05
## Cumulative Proportion 0.997887 0.9991404 0.9996260 0.9999429 1.000e+00
##                           MDS16
## Eigenvalue            4.604e-06
## Proportion Explained  6.561e-06
## Cumulative Proportion 1.000e+00
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CAP1    CAP2     CAP3
## Eigenvalue            0.4027 0.05493 0.009286
## Proportion Explained  0.8625 0.11764 0.019885
## Cumulative Proportion 0.8625 0.98011 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  2.055196 
## 
## 
## Species scores
## 
##       CAP1 CAP2 CAP3 MDS1 MDS2 MDS3
## Dim1                               
## Dim2                               
## Dim3                               
## Dim4                               
## Dim5                               
## Dim6                               
## Dim7                               
## Dim8                               
## Dim9                               
## Dim10                              
## Dim11                              
## Dim12                              
## Dim13                              
## Dim14                              
## Dim15                              
## Dim16                              
## 
## 
## Site scores (weighted sums of species scores)
## 
##         CAP1      CAP2      CAP3      MDS1      MDS2     MDS3
## 1  -0.118345  0.175455 -0.047145 -0.506299 -0.423715  0.09533
## 2  -0.140838  0.003717 -0.302707 -0.387221 -0.320024  0.19578
## 3  -0.103306  0.327975  0.270279 -0.628190 -0.318823 -0.24296
## 4   0.008757 -0.173947  0.898199 -0.589237  0.789919 -0.20068
## 5   0.541369 -0.065379  0.560416 -0.024970  0.049197 -0.32037
## 6   0.543308 -0.004912  0.439088 -0.109438  0.269108 -0.18813
## 7   0.685271 -0.125175  0.720434 -0.172794  0.505845 -0.26832
## 8  -0.368405 -0.840786 -0.581799  0.409126  0.004382 -0.02050
## 9  -0.326275 -0.854491 -0.118938  0.345743  0.065634  0.17900
## 10 -0.447862 -0.569762 -0.479908  0.361270 -0.132091 -0.31818
## 11 -0.271896 -0.855062  0.008803  0.250268  0.264949  0.14661
## 12 -0.363883 -0.842636 -0.568017  0.403435  0.006882 -0.02042
## 13 -0.326080 -0.847913 -0.109894  0.341106  0.062888  0.18602
## 14  0.583087 -0.105603 -0.017661  0.150584 -0.402423  0.38119
## 15  0.552282 -0.099067  0.201860  0.148800 -0.333354  0.02658
## 16  0.682659 -0.146090  0.692094  0.007818 -0.088373  0.36904
## 17 -0.173766  1.053027  0.318616 -0.491137 -0.268441  0.04607
## 18 -0.155959  0.986716  0.683869 -0.490135 -0.299839  0.25781
## 19 -0.184776  1.353849  0.619417 -0.607413 -0.356310 -0.53605
## 20 -0.072462  0.578351  0.794922 -0.392980  0.560962  0.11002
## 21  0.406687  0.358960 -0.891817 -0.118231 -0.441071  1.28428
## 22 -0.526047 -0.062290 -0.245361  0.511732  0.101190  0.06100
## 23 -0.513488 -0.064273 -0.313550  0.508564  0.124313  0.07243
## 24 -0.553768  0.299526  0.216670  0.371002 -0.195805 -0.26911
## 25 -0.482668 -0.375509 -0.390361  0.590368  0.333931  0.25783
## 26  0.614101  0.181112 -0.844084 -0.305500  1.133393  0.20914
## 27  0.194381  0.284271 -0.578667  0.313382 -0.272787 -0.93446
## 28  0.317922  0.429937 -0.934759  0.110348 -0.419536 -0.55896
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##       CAP1     CAP2    CAP3      MDS1      MDS2     MDS3
## 1  -0.2458 -0.44775 -0.1031 -0.506299 -0.423715  0.09533
## 2  -0.2458 -0.44775 -0.1031 -0.387221 -0.320024  0.19578
## 3  -0.2458 -0.44775 -0.1031 -0.628190 -0.318823 -0.24296
## 4  -0.2458 -0.44775 -0.1031 -0.589237  0.789919 -0.20068
## 5   0.5980 -0.09104  0.4327 -0.024970  0.049197 -0.32037
## 6   0.5980 -0.09104  0.4327 -0.109438  0.269108 -0.18813
## 7   0.5980 -0.09104  0.4327 -0.172794  0.505845 -0.26832
## 8  -0.2458 -0.44775 -0.1031  0.409126  0.004382 -0.02050
## 9  -0.2458 -0.44775 -0.1031  0.345743  0.065634  0.17900
## 10 -0.2458 -0.44775 -0.1031  0.361270 -0.132091 -0.31818
## 11 -0.2458 -0.44775 -0.1031  0.250268  0.264949  0.14661
## 12 -0.2458 -0.44775 -0.1031  0.403435  0.006882 -0.02042
## 13 -0.2458 -0.44775 -0.1031  0.341106  0.062888  0.18602
## 14  0.5980 -0.09104  0.4327  0.150584 -0.402423  0.38119
## 15  0.5980 -0.09104  0.4327  0.148800 -0.333354  0.02658
## 16  0.5980 -0.09104  0.4327  0.007818 -0.088373  0.36904
## 17 -0.3329  0.47117  0.2105 -0.491137 -0.268441  0.04607
## 18 -0.3329  0.47117  0.2105 -0.490135 -0.299839  0.25781
## 19 -0.3329  0.47117  0.2105 -0.607413 -0.356310 -0.53605
## 20 -0.3329  0.47117  0.2105 -0.392980  0.560962  0.11002
## 21  0.3833  0.31357 -0.8123 -0.118231 -0.441071  1.28428
## 22 -0.3329  0.47117  0.2105  0.511732  0.101190  0.06100
## 23 -0.3329  0.47117  0.2105  0.508564  0.124313  0.07243
## 24 -0.3329  0.47117  0.2105  0.371002 -0.195805 -0.26911
## 25 -0.3329  0.47117  0.2105  0.590368  0.333931  0.25783
## 26  0.3833  0.31357 -0.8123 -0.305500  1.133393  0.20914
## 27  0.3833  0.31357 -0.8123  0.313382 -0.272787 -0.93446
## 28  0.3833  0.31357 -0.8123  0.110348 -0.419536 -0.55896
## 
## 
## Biplot scores for constraining variables
## 
##                             CAP1    CAP2    CAP3 MDS1 MDS2 MDS3
## isolatewarm               0.2099 -0.9335  0.2908    0    0    0
## temperature28             0.9828  0.1359 -0.1253    0    0    0
## isolatewarm:temperature28 0.8041 -0.1224  0.5818    0    0    0
## 
## 
## Centroids for factor constraints
## 
##                   CAP1     CAP2     CAP3 MDS1 MDS2 MDS3
## isolatecold   -0.09415  0.41864 -0.13043    0    0    0
## isolatewarm    0.07062 -0.31398  0.09782    0    0    0
## temperature6  -0.28450 -0.03934  0.03628    0    0    0
## temperature28  0.51211  0.07081 -0.06531    0    0    0

anova(cap,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for capscale under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = lipids_distBC ~ isolate * temperature)
##          Df SumOfSqs       F Pr(>F)    
## CAP1      1  0.40274 41.1581  0.001 ***
## CAP2      1  0.05493  5.6138  0.034 *  
## CAP3      1  0.00929  0.9489  0.360    
## Residual 24  0.23484                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cap,by="terms",model="direct",perm.max=9999,step=1000)
## Permutation test for capscale under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = lipids_distBC ~ isolate * temperature)
##                     Df SumOfSqs       F Pr(>F)    
## isolate              1  0.06640  6.7859  0.004 ** 
## temperature          1  0.38445 39.2889  0.001 ***
## isolate:temperature  1  0.01611  1.6460  0.189    
## Residual            24  0.23484                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# all functions more or less taken from RDA
# e.g. to get site scores for plotting
#scores(cap,display="sites")[,1:2]